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ABSTRACT 



The immediate attention of the control systems engineer 
is directed to the dynamic behavior of the system under 
study. It is important to study the effects on overall 
system performance of varying one or more parameters (mass, 
inertia, gain, resistance, etc.). It is equally important 
to determine whether a desired dynamic behavior can be 
achieved with any set of values for the parameters — if not, 
redesign is indicated. 

In this thesis a control systems analysis package is 
developed using parameter plane methods. It is an inter- 
active, user-friendly computer aid. Given a characteristic 
equation containing two variable parameters, the output of 
the analysis may be either tabular or graphical, with plots 
of any of the following types: 

1. Constant damping curves as a function of frequency, 

2. Constant frequency curves as a function of damping, 

3. Constant sigma lines (real root lines), 

4. Constant zeta— omega (damping-frequency) curves. 
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I . INTRODUCTION 



The analysis and synthesis of linear feedback control 
systems, or the compensation of same, can be realized by 
three general methods. The first of these can be called the 
integral method. Given a control system, described by a set 
of differential equations, one selects a cost function to be 
minimized with respect to certain variable system parameters. 
The major drawback with this method is the difficulty of 
varying more than one parameter at a time. The second 
method is the Bode frequency response technique whereby the 
system's open loop transfer function is manipulated to obtain 
the desired system response. This method also has its 
inherent weaknesses: difficulty of application to non-unity 

feedback control systems, difficulty in interpreting the 
closed loop transient response in terms of the open loop 
frequency response, and difficulty of varying more than one 
parameter. Third are the algebraic methods. Within this 
category can be included the familiar root locus method. 

Here, a graphical technique is provided by which the set of 
all points which could potentially be made roots are plotted 
in the S-plane. The root locus method is a valuable and 
powerful tool when only one parameter is varied; results 
are less satisfactory for two parameters and of little use 
when three or more parameters are involved. 
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Methods for studying the parameter-root relationship 
when two or more parameters are variable are clearly of 
considerable value. For a linear system, the set of 
differential equations that describe that system can be 
transformed into algebraic equations and manipulated to 
provide a characteristic polynomial. Since the coefficients 
of the characteristic polynomial are deterimined by the 
system parameters, it follows that some relationship exists 
between the value of any parameter and the value of the 
characteristic roots. In reference (1), Mitrovic developed 
an algebraic/graphical method for obtaining the roots of a 
polynomial in terms of two variable parameters. In 
references (2), (3), and (4), Choe, Hyon, and Nutting, 
respectively developed and extended the Mitrovic method to 
the compensation of linear continuous feedback control 
systems. The disadvantage of the Mitrovic method is that the 
variable parameters may appear in no more than two coeffi- 
cients of the characteristic equation, which limits the 
flexibility of the technique. In reference (5), Siljak 
introduced a method for obtaining the roots of a polynomial 
in terms of two variable parameters that may appear in any 
and all of the coefficients of the polynomial. Later, Thaler 
and Towill [Ref. 6] extended this method to the compensation 
of linear continuous feedback control systems. It is from 
the latter work that the ensuing parameter plane equations 
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were developed. General methods of compensation will be 
presented, and an attempt will be made to relate the root 
locus and the parameter plane methods as a set of comple- 
mentary techniques which, when applied in tandem, represent 
the most satisfactory tool to date for designing linear 
feedback control systems. 

The parameter plane method, which works well for two 
variable parameters and which may be extended to three or 
more parameters, is purely algebraic, and the resulting plots 
are valuable aids to analysis. The term parameter plane 
comes from the plot for two parameters — in a rectangular 
coordinate space one parameter will define the abscissa while 
the second parameter defines the ordinate (the S-plane is 
inconvenient for presenting the desired results). Three 
parameters define a 3-dimensional parameter space, etc. For 
design problems it is convenient to think of the algebraic 
calculations as a mapping procedure. By choosing a point on 
the S-plane, the characteristic polynomial acts as a mapping 
function whereby the point may be "mapped" onto the alpha- 
beta plane (alpha and beta are the two variable parameters 
to be used throughout the remainder of this text). The 
relationship between being able to place the roots of a 
polynomial at specific locations in the S-plane and the 
compensation of linear feedback control systems is as 
follows. A feedback control system, including any added 
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compensators which may contain variables, can be reduced to a 
ratio of two polynomials (the closed loop transfer function). 
A specified system response in terms of overshoot, bandwidth, 
settling time, steady-state accuracy, etc., can theoretically 
be obtained by placing a pair of complex conjugate roots of 
the characteristic equation at a specific location in the S- 
plane, while ensuring that the real part of this complex root 
pair (the dominant roots) is smaller in magnitude than the 
real parts of the remaining roots of the characteristic 
equation. The problem of compensation, thus of feedback 
control system design, reduces to one of placing the dominant 
roots of the characteristic equation at the desired location. 
The ability of the parameter plane method to achieve this 
goal will become obvious. 
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II . DERIVATION OF PARAMETER PLANE EQUATIONS 



A linear feedback control system's characteristic equation 
can be expressed as a polynomial of the following form: 



f(s) 



k 

2 a, S - 0, where 
k=0 ^ 



( 2 - 1 ) 



( k=0 , 1 , . . . ,m) are real coefficients 
S = -a+ju) = -|o)+ jo) /l-^^ 

is the undamped natural frequency and 
is the relative damping coefficient 
In reference (5) it is noted that S may be represented by 
the following: 






( 2 - 2 ) 



where 

T^(-U = (-l)^Tj^(0 and Uj^( - 4 ) = ( -1 )^^^Uj^( O • 

Tj^(^) and Uj^(^) are Chebishev functions of the first and 
second kind respectively. Values of zeta and omega will be 
considered such that 0<^<! and 0<o)<oo. Values of T^ and Uj^ 
are tabulated in .various appendixes. More important to digital 
computer analysis, they can be obtained from the following 
recursive relations: 
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0 



- 



- "V«> ^ = ° 



(2-3) 



Here, Tq(4)=1, T^(4) = C, Uq(O= 0, U^(0 = 1. Substituting 
equation (2-2) into (2-1) and setting the real and imaginary 
parts to zero independently, one obtains: 



m 



1 a co'^T (-0=0 
k=0 ^ 



(2-4) 



m 



2 a (-0=0 
k=0 ^ ^ 



Employing equations (2-3), one obtaihs from equation (2-4): 



m 

2 

k=0 






0 



k k 

Z (-l)^a. (O = 0 
k=0 ^ ^ 



(2-5) 



Now consider the coefficients a, of the characteristic 

k 

equation (2-1) as linear functions of the variable system 
parameters, a and (3, as follows: 

^k " ^k“ ^k^ ^ ^k 
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Using this relation for a^, equations ( 2-5 ) become; 

= 0 

^ °2 = ° 

where 



(2-7) 



k k 

B = S (-1)"^ b w'' U 
-L k=0 ^ ^ 



m u V 

B„ = E (-1) b,co U, 
^ k=0 ^ ^ 



= 



m 

Z 

k=0 



, - . k , k tt 

(- 1 ) 



C„ = 



ra 

Z 

k=0 






(2-8) 



D. 



m 

= Z 
k=0 



v" Vi 



D„ = 



m 

Z 

k=0 



. T X k , ,k 

(-1) d^o. 



Since equations ( 2-7 ) are linear in the two unknowns alpha and 
beta, Cramer's rule may be applied to obtain: 



Equations ( 2-9 ) are now functions of zeta and omega. Hence, 
by fixing either zeta or omega and varying the remaining 
parameter, the constant omega or constant zeta S-plane 
contours respectively can be mapped into the real domain 
of the alpha-beta or parameter plane. 

In reference (5) the following relationships are noted: 







0 




0 



(2-10) 





and are related to the Chebishev functions by: 






( 2 - 11 ) 
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By employing equations (2-10) and (2-11), one obtains 
(proceeding as before); 



m 

S 

k=0 






0 



m 

2 

k=0 






0 



( 2 - 12 ) 



Combining equations (2-6) and (2-12) with Cramer's rule, one 
again arrives at equations (2-9) where the following expressions 
now apply: 



B 



1 



m 

2 

k=0 



\ ‘^k-l 



0 



B. 



m 

2 

k=0 



b,, Q, 



0 



c 



1 



m 

2 

k=0 



=k “^k-i 



0 



m 

2 

k=0 



Q, = 0 



(2-13) 



D 



1 



m 

2 

k=0 



^k \- 



0 



D. 



m 

2 

k=0 



Q, 



0 



Equations (2-9) and (2-13) are useful for mapping constant 
zeta-omega curves from the S-plane into the parameter plane. 

As will be demonstrated later, these curves play an important 
role in dominance considerations. 

If the complex variable S is substituted in equation (2-1) 
by letting S =— a, where sigma corresponds to values of S along 
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the real axis, then according to equation (2-6) the 
characteristic equation (2-1) becomes: 

^ kk ^ kk^ kk 

0! 2 (-l)%,a^ + ;3 2 (-1) c,CT^ + I (-l)^d = 0 

k=0 ^ k=0 ^ k=0 ^ 

(2-14) 

The above expression represents a straight line in the alpha- 

beta plane for a given value of sigma. Hence a point on the 

real axis in the S-plane maps into a straight line in the 

alpha-beta plane. In addition, for given values of alpha, 

beta, and sigma which satisfy equation (2-14), the 

characteristic equation (2-1) must have a real root at minus 

sigma. On the constant zeta and omega curves previously 

defined, for certain values of alpha and beta (say, for 

values obtained from equations (2-9) for given values of 

zeta and omega) the characteristic equation will have a pair 

2 

of complex conjugate roots at S = -§cj + jco/ 1-5 

The significance of the above discussion is that by 

applying equations (2-9) and (2-14) one can, for a specified 

value of zeta, omega, and sigma, compute the value of alpha 

and beta such that the characteristic equation will have a 

2 

pair of roots at S = ^”^1 * 

remaining roots of the characteristic equation can then be 
calculated by dividing out the two known or specified 
roots. This method, where zeta, omega and sigma, or simply 
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zeta and omega are specified, and where the computations for 
alpha and beta are done algebraically, will be referred to as 
the algebraic parameter plane solution. 

To solve the problem in general for all values of zeta, 
omega, and sigma, it becomes necessary to plot a family of 
parameter plane curves for various values of zeta, omega, 
sigma, and if desired, zeta-omega. On the resulting 
parameter plane plot one can, by choosing an operating point, 
graphically read from the curves the values of alpha and beta 
and, hence, the values of the m roots of the corresponding 
m^*^ order characteristic equation. This latter method will 
be referred to as the graphical parameter plane solution. 

The algebraic solution has the advantage that the labor 
of plotting the curves can be avoided, but the disadvantage 
remains that without the curves it is difficult to pick the 
optimum values of zeta and omega so as to ensure dominance 
while still meeting the system specifications. The graphical 
solution has the advantage that one has a "picture" of the 
way the characteristic roots move about in the S-plane as 
alpha and beta are varied. This enables one to choose the 
values of alpha and beta corresponding to the best values of 
zeta, omega, sigma, and zeta-omega for all roots of the 
characteristic equation. This feature of the parameter 
plane points out a strong justification for attempting to 
obtain the parameter plane curves. And with the employment 
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of a digital computer and an appropriate algorithm to realize 
the parameter plane curves, the advantage of the algebraic 
method becomes muted. 

Recursion methods (equation (2-3)) are by no means the 
only methods of producing algebraic and graphical parameter 
plane data. Thaler and Karmarkar [Ref. 7] describe a matrix 
solution to the parameter plane problem. Essentially, a 
matrix of coefficients may be manipulated to obtain the 
following general form: 



— 






















^1 


2 

— U) 


0 


0 


. . 0 


<^1 


®1 




a 


^2 


^2 


-2^w 


2 

— w 


0 


. . 0 


^2 


®2 




■ 6 


• 


• 


-1 


-2^« 


2 

-U) 


. . 0 


• 


• 




m-2 


• 


• 


• 


-1 


-2?w 


. . 0 


• 


• 




• 


• 


• 


• 


• 


• 


2 


*^m-2 


• 




m-2 

^1 


• 


• 


• 


• 


• 




H 2 

a -1 
m-1 


• 




1 


^m 


c 

m 


0 


0 


0 


. . -1 


d 

m 


e 

m 




as 
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where bj^, , d^, w , and are as described before, and: 

e, (k=l, m) are the coefficients associated with the 

non-linear alpha-beta product terms 

R. is the sum of the m-2 roots of the polynomial character- 
istic equation taken one at a time 

m-2 

R „ is the sum of the m-2 roots taken m-2 at a time 
m— 2 

Further, by the application of appropriate row operations, 
this matrix may be reduced to the following row-echelon form: 



1 


0 


0 






0 


k, , 


k-, „ 




a 


0 


1 


0 


• 


• 


0 


11 

^21 


12 

^22 




6 


0 


0 


1 


• 


• 


0 


^31 


^32 




m-2 

m-2 


• 


• 




• 


• 


«• 


• 


• 




' 


• 


• 


• 


• 


• 


• 


• 


• 




♦ 


• 


• 


• 


• 


• 


• 


• 


• 




1 


0 


0 


0 






1 


k . 


k „ 




aB 














ml 


m2 ' 
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For the case when all 
equation are linear, 

( k=l , . . . ,m ) = 0 , and 



coefficients of the characteristic 
i.e. k=l , . . . ,m) = 0, then 



^ = -K 



11 



6 = -K 



21 



Us obtained from the first two rows of the matrix equation. 

One should note that in arriving at equations ( 2-15 ) , 

2 

approximately m row operations are required for the row- 
echelon matrix formulation for each point of the parameter 
plane curves (e.g., each time either zeta or sigma are 
varied). Compare this with the approximate m calculations 
required to obtain the recursion equations of the previous 
chapter, and the matrix method becomes relatively inefficient 
for larger order systems. 

One should not, however, discard the matrix approach 
entirely. For small order characteristic equations, this 
technique compares favorably with the recursion method. 

And when the variable parameters are non-linear — when one 
must deal with alpha-beta product terms — the matrix approach 
affords a more direct method of obtaining the alpha-beta 
pairs. ^ Whether the recursion or matrix method is utilized 



'From equations (2-15), one obtains the two quadratic forms 



S 2 ® ■^(''21 ‘^ 12-'^11 









from which alpha and beta are easily derived. 
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should depend on the inclusion of alpha-beta product terms; 
ultimately, it is a matter of personal preference. 
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III. APPLICATION OF THE PARAMETER PLANE METHOD 



A. ALGEBRAIC SOLUTION 

In this section it will be assumed that the system 
performance specifications can be met by placing a pair of 
complex conjugate roots at a specific location (i.e., by 
choosing appropriate values for zeta and omega). If after 
computation of the necessary values of alpha and beta to 
locate the roots as desired it is found that these specified 
roots are not dominant , then either a different value of zeta 
and/or omega must be used (possibly at the sacrifice of some 
performance measure), or a different method of compensation 
will have to be attempted. In a later section a method will 
be addressed whereby the dominancy requirement may be^ 
achieved. 

1. Feedback Compensation 

For a unity feedback control system, let 



G = 



K 

e(S) 



K 

s"*+e . .+e, 

m-1 L 



(3-1) 



where K is the forward path gain (a variable) and e(S) is a 
polynomial in S representing the poles of the open loop 
transfer function of the uncompensated system. In equation 
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(3-1), L corresponds to the system type--for a type 0 
system, L=0, for type 1, L=l, etc. The system's error 
coefficient is defined as: 

K = lim (3-2) 

® S-^’O 

where is t'he open loop transfer function of the 

compensated system. 

a. Tachometer Plus Acceleration Feedback 

In order to achieve the system performance 
specifications, a feedback compensator must be introduced. 
Let 

H = K^S + K 
t a 

The resulting compensated system's characteristic equation 
becomes : 

e(S) + K(K,S+K S^) = 0 (3-3) 

U cl 

and by expanding e(S), equation (3-3) becomes: 

s”^+e . . + (e„+KK )S^+( e, +KK^ )S+e +K = 0 (3-4) 

m-1 2 a lt o 

where L is zero for a type 0 system (the most general case). 
The following results also apply to a type 1 system if e^ is 
set to zero, and, similarly, for a type 2 system if both 
and e^ are set to zero, etc. Combining equations (3-2), 
(3-3), and (3-4) the error coefficient becomes: 
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K 



(3-5) 



K 

e 



lim 

S-^0 



S°K 



e(S)+K(K.S+K S ) 
t s. 



e +K 
o 



or for a type 1 uncompensated system: 
K = K 

or for a type 2 uncompensated system; 

K - K 

e - KK^ 



(3-6) 



(3-7) 



Note that if the uncompensated system is type 2, the 
compensated system would be type 1 if tachometer feedback 
or tachometer plus acceleration feedback is used. 

In the compensated system's characteristic 
equation (3-4) let alpha = KK and beta = KK . Equation 
(3-4) then becomes: 

+ e T ^ +. . . + (er,+oi)s^+(e +6)S+e. +K = 0 
m-i z 1 o 

Recalling equation (2-6) where in general the coefficients 
of the characteristic equation are of the form: 



a 



k 



bka+c^6+d^ 



and letting m=k, then from equations (2-8) one obtains: 
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^2 = " 



= -cjU = 0 
1 o 



C^ I — — wU^ — —01 



(3-8) 



^ k k ^ k k 

D = E (-1) d <0 U D = E (-1) d 

1 K K i ^ k k 



since U =0 and U^=l. From equations (2-9) one derives: 
o 1 

k 

C D -C D ^ ( ~1) ‘^k^k-l m 

a = 12 ^2^1 ^ _k=0 ^ ^ k ,k-2 

B,C_-B_C, 3 I-* '^k ^k-1 

12 2 1 -0) k=0 



(3-9) 



m b- b- 1 

k=0 



If alpha and beta are linear functions of K, the forward 
path gain, one can use the steady -state error specification 
to define K in terms of alpha and/or beta. Since zeta and 
omega were assumed to be specified, then from equations 
(3-9) one may solve for alpha and beta. From this, and 
are readily determined. 

Example 3-1 

The system of Figure (3-1) is to be compensated by using 
tachometer plus acceleration feedback . The system 
specifications are as follows: 
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/K 



K 

5(S+1)(S+2)~ 




Figure 3-1 



1. Complex roots corresponding to C=0.7 and “*=10. 



2. K >6 
e— 

From equation 



e 2+KK^ 



(3-2): 
> 6 



From which K_>12 + 6KK^. The compensated characteristic 
equation is 



S^+S^(3+KK )+S(2+KK^)+K = 0 
a t 



(3-10) 
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Letting alpha = KK and beta = KK , equation (3-10) becomes; 

8L L 

S^+S^(3+cv)+S(2 + /3)+K = 0 (3-lOa) 

From equations (3-8): 

= 100 = 140 

Cl = 0 Cg = -10 

D = -1100-K = -1120 

and from equations (3-9): 

^ _ lO(-llOO-K) 140(-1100-K)+56000 n 

-1000 ’ ^ -1000 ^ ^ 



From the steady-state accuracy specifications, then, it is 

necessary that K>12+6/3; let K=12+6/3. From equations (3-11) 

it is found that /3 = 623, hence K=3750. Therefore, Q!=48.5, 

and since 0i=KK and /3=KK, : 
a t 



K = 

a 3750 



0.0129 



623 
t 



3750 



0.1661 



The compensated system's characteristic equation becomes 



S^+51.5S^+625S+3750=0 



(3-12) 




Dividing equation (3-12) by this quadratic, the remainder 
is S+37.5. Since zeta'omega of the desired roots = 7<<37.5, 
the complex roots are dominant and the problem is solved, 
b. Tachometer Feedback Only 

Let H = K^S. The characteristic equation of the 
compensated system becomes: 



. .+e2S^ + (e^+a)S+eQ+3=0 



Proceeding as in the previous example, one obtains: 



B 



1 



0 



Br, - -0) 



C 



1 



-1 



C 



2 



0 



m 




m 
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and 



a 




B = S 
k=0 



m 




(3-13) 



For a specified value of zeta and omega, alpha and beta can 
be obtained from equations (3-13). The error coefficient is 
then determined directly from equations (3-5), (3-6), or 
(3-7). Thus, the error coefficient is fixed for a given 
value of zeta and omega, and if this parameter is to be met, 
the values of zeta and omega may require adjustment. One 
possible approach might be to fix zeta at some value, 
whereby from the given and equations (3-5), (3-6), or 
(3-7) alpha could be computed. Equations (3-13) could then 
be solved for, first, omega and then beta. The calculations 
would prove tedious, however. 
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Example 3-2 



The same system as used in example (3-1) will be studied 
here, this time with tachometer feedback alone. The same 
system performance specifications are to be met, namely, 

K^> 6 , zeta = 0.7, and omega = 10. The compensated system’s 
characteristic equation becomes: 

S^+3S^+(2+KK^ )S+K = 0 (3-14) 

Letting Q!=KK^ and 13 = K here, equation (3-14) becomes: 

S^+3S^(2+Q!)S+/3 = 0 

From equation (3-13) it is found that: 
a = -2+30(1. 4)-100(0. 96) = -56 

Since alpha is negative, it is seen that positive tachometer 

feedback is required. Further, it is found that the 

remaining root (when equation (3-14) is divided by 
2 

S +14S+100) is positive; the system is unstable. Hence 

the desired system specifications cannot be met with 

tachometer feedback alone. 

c. Acceleration Feedback Only 

2 

Let H = K S . The characteristic equation of 
a 

the compensated system becomes: 
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s'^+e ^+...+(e„+KK )S^+e,S+e +K = 0 

m- 1 2 a 1 o 



Proceeding as before, where now Q!=KK and /3=K; 

a 



Bf = = 



^2 = “ ^2 



C, = 



= -1 



c„ = 0 



m , , 

D = S (-1)X<^^U 
■L k=0 ^ ^ ^ 



m 



D = 2 

k=0 ^ 



Solving for alpha and beta yields: 



a = 






1 f , ,.k , k-2., 

^ ^ k k 



m 

/? = 2 
k=0 



. , N k , k„ 

(-1) dj^o* U 



k-1 




m 

2 

k=0 



. T ^k , k„ 

(-1) d^- 



(3-15) 



Calculations for alpha, beta, and are performed in the 
same manner as with the preceding tachometer feedback 
example . 

Example 3-3 

The same system of examples (3-1) and (3-2) will now be 
compensated using acceleration feedback alone. As before, 

K >6 , zeta=0.7, and omega=10. Therefore K ~ and the error 

C ' 0 ^ 
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coefficient is unaffected by the acceleration feedback. Hence 
one can conveniently choose K=12 to meet the specifications. 
The compensated system's characteristic equation becomes: 

S^+(3+KK )S^+2S+K = 0 (3-16) 

a 

If K in equation (3-16) is set equal to 12 as prescribed, only 
one parameter remains and the parameter plane equations 
produce an indeterminate solution. If K is left as the 
variable beta, then equation (3-16) becomes (after the usual 
substitut ions ) : 

S^+(3+Qr)S^+2S+/3 = 0 

By employing equations (3-15), one obtains o;=4 and /3= -700. 
Since beta is negative, it is concluded that the desired 
roots (i.e., desired values of zeta and omega) cannot be 
realized using acceleration feedback alone, and of course 
neither can the desired error specification be obtained. 

One would therefore choose an alternate method of 
compensation . 

If one chooses to use feedback compensation then 
perhaps tachometer plus acceleration feedback might be 
attempted first using equations ( 3-9 ) and the appropriate 
steady-state error specification. If the specifications cannot 
be met in this manner, then it follows that neither 
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t ac lK>inot or nor acci- 1 (.,■ r at i on foocibai.'k alotu- will .-,u t f i (.m- . In 
this case either the system's spt.'c i f i cat i ons must be eased t;r 
another type of compensation must be utili-^'.ed. If it is 
found that the specifications are achievable with the combined 
tachometer and acceleration feedback, then, if desired, 
equations (3-13) and (3-15) can be employed lo investigate 
the feasibility of tachometer or a<.:ce 1 era t i en feedback alone, 
respect i ve ly . 

d. Case For Which Feedback Is Not Available Near 
The Forward Path Amplifier 






K 


•^(T)-> 




Y- 



S(S + 1)(S+2) 




Figure 3-2 

Figure 3-2 shows a system similar to that used 
in example (3-1) except that now the feedback is inserted at 
the output terminals of the amplifier represented by gain K. 
This illustrates a system for which it may not bo possible 
or practical to access the input terminals of the error 
detector. This problem will be solved by means of an example. 
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Example 3-4 



As before, the same system specifications are to be met, 

» zeta=0.7, and omega=10. The characteristic 
equation becomes: 

S^+(3+K )S^+(2+K^ )S+K = 0 
a t 

Letting o;=K and j3=K : 

3 . u 

S^+(3+a)S^ + (2+/3)S+K = 0 (3-17) 

Comparison of equations (3-17) and (3-lOa) show that they 
are identical, that is, the solution obtained for alpha and 
beta in example (3-1) applies. There, alpha and beta were 
found to be 48.5 and 623, respectively, while K=3750. For 
the present example no further computations are necessary 
to find K and K , since they are now the parameters alpha 

3. L 

and beta. This points out an important advantage of the 
parameter plane method, namely, that the solutions depend only 
on the characteristic equation and not on the system from 
which the characteristic equation was formed. This principle 
can similarly be applied to control problems involving 
tachometer or acceleration feedback alone. 

2 . Cascade Compensation 

For a unity feedback control system let G have the 
form of equation (3-1): 
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K 



( 3 - 1 ) 



m-1 



+...+ 6^3 



L 



where again K is the forward path gain (a variable) and 
e(s) is a polynomial in S representing the poles of the 
open loop transfer function of the uncompensated system. 
The letter L again indicates the system type. If in order 
to satisfy the system's requirements a cascade compensator 
is required, then let; 



With a d.c. gain of unity, placement of this compensator in 

the forward path will not affect steady-state accuracy. 

With Gq as indicated here, the values of P and Z are 

computed to obtain the desired system response. If P is 

less than Z, a lag network is required and the factor of 
p 

2 of the compensator is inherently present due to the 
physical nature of the compensator (usually an R-C network). 
In this case all forward path amplifier gains can remain 
unaltered to meet the specified accuracy demands. If, 
however, the computed value of P is greater than Z, a 
lead network is required and the compensated system's 

p 

forward path gain must be raised by the factor of ^ 'to 
meet the accuracy specifications. As the physical nature 

p 

of the lead network is such that the factor 77 is not 



P(S+Z) 
C Z(S+P) 
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inherently present, this factor must be provided either by 
adding an amplifier in cascade with the lead network or by 
raising the gain K of the existing amplifier as required to 
achieve steady-state accuracy. 

Continuing, the compensated system's forward path 
transfer function is: 

P _ P P _ K P S+Z _ K 

^CC e(s) ' Z * S+P S+P * e(s) 



Applying the definition of the error coefficient one 
obtains : 



K 

e 



lim S^‘ 
s->-0 



K 



e(s ) 



y(s+y) 

(S+P) 






and again assuming a type 0 system where L=0 , the compensated 
system's characteristic equation becomes: 



e(sKS+P)+KY(S+y) = 0 



or after expansion: 



i t / I \ I / L \ ^ I 

S +(P+e^ ^ )S +(Pe^ +e )S +... 

m— 1 m-1 m— 2 



+(Pe„+e, )S +(KY+Pe,+e )S+P(e +K) = 0 
2 1 1 o o 



(3-18) 
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Letting a=P and S=Y, equation (3-18) becomes: 




+ ( e2^‘*'e^ )S^ + ( K3 +e^ot+eQ ) S+ci ( Bq+K ) = 0 



(3-19) 



Comparison of equation (3-19) with the general form of the 
characteristic equation as specified in equations (2-1) and 
(2-6), it is apparent that K=m+1, bQ=eQ+K, CQ=dg=0 , t>^=e^, 



variable beta represents the pole-to-zero ratio of the 
cascade compensator. The S-plane can be divided into regions 
where lag compensation or lead compensation is needed. By 
mapping of variables in the above manner, the parameter 
plane can effectively be divided into corresponding regions 
above and below the straight line 6=1. Then, for values of 
beta less than one a lag network is required and for beta 
greater than one a lead network is needed. In addition, if 
beta is less than 0.1 or greater than 10, a multiple lag or 
multiple lead network, respectively, is required. 



Cq^ K, ^0^ ^2 ^2^ ^2 ^ ^ ^2 ^1^ etc. 

It is important to note that the parameter plane 



Based on equations (3-19) and (2-8) it is found that: 



B^ = -( eQ+K) + to^e2+ •••■•■( -1 ) 



k-2, k-2 

GO 






c 



1 



0 
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D ^ = 01 ‘ 



1 , . / 1 sk-2 k-2 

e^ + . . . + (-1) 

1 m-z 



^k-3 ^ ®m-l '^k-2 



. f , sk k,, 
+ (- 1 ) - 



®2 ^ -‘“6^ + " 



^2^2^ 






C 



2 



-0)K 



D_ = -o»e +" 
2 o 






, . , sk-2 k-2„ 

.+(-1) e^_2- U^.2 



+<-l)‘''^e 



m-l 



>-lUj^_l + (-l)^“>^Uk 



(3-20) 



and from equations (2-9) it is found that: 



a = 




P = 






(3-21) 



For a type 1 system, e^ in equations (3-20) is set equal to 
zero, for a type 2 system e^=e^=0, etc. On the basis of 
equations (3-20) and (3-21) a cascade compensator can be 
designed. 

Example 3-5 

For the system of Figure (3-3) it is desired to design a 
cascade compensator which places a pair of roots at zeta=0.5 
and omega=l. The error coefficient should be 50. 
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Figure 3-3 

It is apparent from Figure (3-3) that 
K=500. The characteristic equation is: 

S‘^+(8 + P)S^+( 17 + 8P)S^ + (10 + 17P+Ky)S + P( 10 + K)=0 

and by the substitutions of -i=P and p = i: 

S^ + (8 + 0')S^+( l7 + 8ii)S“ + ( l0 + l7a + 500d)S + 5l0a' = 0 

Applying equations (3-20) and (3-21) one obtains: 

= -503 = -9 

= 0 = -500 

= 9 ^2 = 6 

= 0.0179 = P ^ = 0.0117 = y 
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p 

and since Z=1.529. This is a lag network for which the 

factor 0.0117 is inherent in the R-C filter design. 

Although treatment will not be presented here, the 
algebraic application of the parameter plane technique can 
be readily applied to combination cascade and feedback 
compensat ion . 

B. DOMINANCY OF THE SPECIFIED ROOTS 

In the preceding section nothing was done in the 
calculations to make the specified roots a dominant pair. 

As mentioned earlier, the ability to predict a system's 
response on the basis of the location of a pair of complex 
conjugate roots was based on the assumption that the 
magnitude of the real part of the specified roots was much 
less than that of the remaining roots of the characteristic 
equation. In practice, if the real part of the specified or 
primary roots is one half to one fifth or less of the real 
parts of all secondary roots, the system is said to be 
dominant in the primary roots. In many cases the system 
will still meet the specifications even if two pairs of complex 
roots have the same real part, provided the zetas for both 
pairs of roots meet the specifications, and the undamped 
natural frequencies are such that the component time responses 
are not highly additive. Further, even if there exists a 
characteristic root whose real part is closer to the origin 
than that of the primary pair, the presence of a closed-loop 
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zero could make the residue of the close-in root negligible 
as compared to the primary root. If possible, however, one 
usually attempts to make the real parts of all secondary 
roots 1 ai'ge in magnitude. 

In the preceding examples it should be pointed out that 
in many cases there were actually three and possibly four 
variable parameters. For instance, the forward path gain was 
usually a fixed value in the computations so as to meet the 
minimum steady-state accuracy requirements. There is, 
however, no reason why the gain cannot be raised above the 
minimum value, thus permitting a third degree of freedom. 

When cascade and feedback compensation are used simultaneously, 
the forward path gain and tachometer gain become the third 
and .fourth parameters. 

Recall that the system characteristic equation has the 

k 

form f(S) = 2 aj^S =0, where realize the 

k=0 

system specifications, one places a pair of complex roots 
at S=- ^ which implies that S^+2C^aj^S+oo^=0. 

Since and are known, the quadratic can be divided 
out of the characteristic equation, leaving a polynomial 
which contains the secondary roots of the characteristic 
equation. Since only two of the variable parameters were 
used in fixing the primary roots, the remaining parameters 
will appear in the coefficients of the quotient polynomial, 
and it is these parameters that can be varied to achieve 
dominance . 
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Instead of division to obtain the quotient polynomial, 
coefficients of like powers will be equated to achieve a 
system of equations. Let the quotient polynomial be 
given by: 

" k 

f,(s) = I f.S^ = 0 (3-22) 

^ k=0 ^ 

where n=m-2 , i.e., equation (3-22) is of order two less than 
the characteristic equation. Applying equations (2-1), 
(3-22), and the quadratic it follows that: 



2 _ 2 ^ k k 

(S^+2C CO s+w^)( I f,S^) = I a.S^ 

-L -L k=o ^ k=0 ^ 



(3-23) 



Taking a^^=l and equating coefficients of like power: 



a = f =1 
m n 






m-2 n-2 ^1 1 n-1 1 



^2 = 



a = 2 ^ CO f +f jj 
1 1 1 o 1 1 



(3-24) 



^0 = ^o“ 1 
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Equations (3-23) and (3-24) can be solved for the coefficients 
f in terms of the coefficients a. The results will be 
applied to the following cases: 

Case of k=3, n=l 
Equation (3-23) becomes: 



(S^+2 5^a)^s+03^)(f^S+fQ)=S^+a2S^+a^S+aQ 



Equating coefficients of like power one obtains: 



a 



3 



f 



1 



1 






Solving for the coefficients f results in: 



f 



1 



1 




(3-25) 
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Case of k=4 , n=2 



Proceeding as before, the a coefficients become: 



= f2 = 1 



ag - 






= 






^0 



f <0? 
o 1 



(3-26) 



When solved for the coefficients f, equations (3-26) yield: 



^2 = 1 



CO, 



f. = 



0 



^3-2^i"i = 


2C 


^0 




a CO 

-31,,^ 


co2 ■ 


2«l“l 





0 



2^lC0i3 



1 , ^^ 1 ^ 0 . 
CO “ ^ ‘^l 






(3-27) 
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Case of k=5 , n=3 



Similarly : 



^5 ^3 ^ 



»4 - 2 5i"i+t2 



^3 =“l 



^2 =“l V2«l"lV^0 






“l*0 



and : 



fo = 1 



f., = 



2 

a„ 4§ 2^ a 

^4-2 4 ^1 = ~2 ^ 



^0 “l 



f. = 



^^1 1 ^^1 1 1 

^-1 2? a^ 2 2 

Cj— —3 ^ ^3 ~ ^^l“l^4 '^^l^l “ “ 



f = ^ 
0 0)2 



43 



Although the coef f icients o f f have been derived for only up to 
the fifth order case, they can easily be obtained for higher 
order cases if necessary. 

Example 3-6 (Third Order Characteristic Equation) 






(S+P/ ) 




K 


S+P 




S(S+4) 



-> 



Figure 3-4 

Design a cascade compensator for the syst-em of Figure (3-4) 
to obtain: 

1. Characteristic roots at zeta=0.5 and omega=40. 

2. K >250. 

e— 

3. The specified roots are to be dominant. 

The characteristic equation of Figure (3-4) is: 

+ (4+P)S^ + (4P+Ky)S + KP = 0 
or 

+ (4+«)S^ + (4«+Ki3)S +K« = 0 
where “=P and /? = y. 
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Here G s^+4s’ 

Application of equations (3-20) yields: 




= -K+O) = -K+1600 



B_ = -4&) + CO U_ = 1440 



= 0 = -coK = -40K 

= 40)^ - (o^U^ = -57600 = ^co^U^ - co^U^ = 6400 

From equat ions ( 3-21 ) are obtained; 

^ _ 57600 
-K+1600 

/? _ 1440(-57600)-(-K+1600)(6400) 

_40K( -K+1600) 

For any value of K, equations (3-28) will produce values of 
alpha and beta that provide characteristic roots at zeta=0.5 
and omega=40. However, only certain ranges of K will meet 
the steady-state error and dominance requirements . To satisfy 
the error consideration it is necessary that K>1000. Since 
f^(s) is of order one, equations (3-25) apply and; 



f 

f 



1 

0 



1 




( ) 



- 
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From the characteristic equation it is seen that 



^2 = 

a^ = 4 a + K 6 



aQ = Ka 



The real part of the specified roots is |co=20. Arbitrarily 
choosing a dominance factor of five, the dominance criterion 
becomes: fg > 5|oo = 100. To satisfy this requirement, the 



simplest form of fg will be chosen, namely 



then seen that f = 



Ka 



57600K 



OJ 



1 36K 



:> 100 . 



o 1600 1600C-K+1600) 1600-K 

This requires that K>1176.5. Since KM176.5 also satisfies 

the error specification, a value of K=1180 is arbitrarily 

chosen. Using this value of K, it is found that: 

a=l37, 8=4.3 and f^=101. As a check, from the expression 



fQ = (4+137)-2(0.5)(40) =• 101 
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Example 3-7 (Fourth Order Chan 



i c E(ju a L i cjn ) 




K 



A 



^ (S+0.5)(S+1)(S + 5)(S-M0) 



2 I 

KqS + K^Sj^ 



1 



J 



Figure 3-5 



Compensate the system of Figure (3-5) using tachometer plus 
acceleration feedback to obtain: 

1. Characteristic roots at zeta = 0.5 and omega = 2. 

2. K >12. 

e— 

3. Dominance of specified roots. 

The characteristic equation of Figure (3-5) is: 



S^+16.5S^+(73+a)S^+(82.5+6)S+25+K = 0 



where u=KK 



a 



and S=KK^. 




K 



+ 16. 58^+7 33“ + 82 .53 + 25 
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By inspection 6^=25, e^=82.5, 02=73, 02=16.5, e^=l, 62=!, 

62=0, and U^=-l . When equations (3-9) are employed one finds: 

^ _ K-135 a _ K-24 

4 ’ 2 

From the characteristic equation it is seen that a^=l, a2=16.5 

a2=73+“, a^=82.5+6, aQ=25+K. Since the quotient polynomial 

2 

f^(s) is a quadratic, i.e., S +f^S+f2=0, equations (3-27) apply 
It would be desirable that, from a dominancy viewpoint, 
f^>5 = 5. However, looking at the dominancy equations 
for this case (equations (3-27)), it is seen that one of 
the several expressions for f^ is f^=a2-2|^“^, which is a 
fixed constant even though the remaining expressions for f^ 
involve one or more variables. Thus, f^=14.5. Noting the 
most simple expression for fp. in equations (3-27), f_= — 

Now, since f^=14.5>5, a dominant situation already exists. 
However, the system's performance can be further improved 
by choosing appropriate values of zeta and omega for the 
secondary roots. From the error specification it is 
necessary that -^>12 or K>300. Now f ^ ( s )=S^+14 . 5S+— ^ 
or f^(s)=S^ + 14.5S+6.25+0.25K. For K=300 , f^(s)= ^ 

S^ + 14 . 5S+81. 25 . Therefore, 212^2“^“^*^’ ^"^=^ 1.25 implying 
Then zeta=0.806. These appear to be reasonable 
values for ^2 ^2 secondary roots taken aJone 

would produce less overshoot and a smaller settling time 
than the primary roots. Using this value of K one obtains 
for alpha and beta: 
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a = 41.25 



6 = 138 



and since a=KK, and 6=KK , K is found to be 0.1375 while 

L 3. L 

K =0.46. 
a 

As an added bonus of this method, all the roots of the 
characteristic equation are now known and the system's 
time response could be computed if desired. 
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IV. PROGRAM DESCRIPTION 



The goal in developing any computer aided design program 
should be twofold; (1) provide the user with a single, 
easily understandable, easily usable and comprehensive 
engineering tool, and (2) dramatize its efficiency above 
that of other currently available methods. Through the examples 
of the following chapter the second goal will be demonstrated. 

It is first desirable to reveal the methodology and internal 
structure of the parameter plane curve program — as a 
consequence it is hoped that the first goal will be affirmed. 

A. THE PROGRAM 

The parameter plane curve — generating program, or 
"program" as it will be called henceforth, consists of a 
large driving routine which includes all necessary 
calculations with which to generate the curve data, and 
several supporting subroutines (i.e., curve plotting, root 
solving, data saving, etc.). This entire package is 
included as a user-selected option within another controls 
system computer aided design package. Among other options, 
the latter CAD program includes a root-locus analysis-- 
as mentioned earlier, the usefulness of either the parameter 
plane or root locus technique for design of a controls 
system is somewhat limited, but in combination their 
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effectiveness is synergistic (Chapter V will assert the dual 
roles of the root-locus and parameter plane methods). 
Facilities available within the program are many; the 
major options are: 

1. Plotting of constant zeta curves, with zeta as a 
function of omega. 

2. Plotting of constant omega curves, with omega as a 
function of zeta. 

3. Plotting of constant sigma (real root) curves. 

4. Plotting of constant zeta-omega curves. 

5. Tabular output. 

6. Rescaling of the plots. 

Input of certain data is required to enable the program. 
These inputs include: 

1. Starting value of co . 

n 

2. Decades of w to be considered. 

n 

3. Number (and values) of constant zeta, omega, sigma, 
and zeta-omega curves. 

4. Coefficients associated with the constant, alpha, and 
beta terms of the characteristic equation. 

Each of the basic program option areas will be described in 

appropriate detail, as well as their interaction with the 

input data. 

1 . Constant Zeta Contours 

In practice design specifications for control 
systems are given in terms of percent overshoot, settling 
time, error constraints, etc. A value of zeta can be 
associated with the first of these specifications — that is. 
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a given percent overshoot requirement can be related to a 
specific value of zeta. Given a specific zeta value, the 
program calculates the alpha and beta coefficients of the 
characteristic equation by holding the zeta value constant 
and varying the value of omega. The limits within which 
omega is varied are defined by the user's choice of the 
initial value, and the number of decades of omega to 
be considered. 

From the nature of the mapping process, it is clear 
that when the contour of the coefficient plane passes through 
a designated point (M-point), the original mapping contour 
on the S-plane passes through a point which is a root of the 
characteristic equation. The zeta value chosen for the 
contour is then the zeta for the root. The value of omega 
associated with the M-point is' the radial distance from the 
origin of the S-plane to the root. Thus, a complex root 
is determined when the M-point lies on a constant zeta curve 
of the parameter plane. The value of this root and its 
complex conjugate is: 

S = -^co+jo)/ 1-1^ 

If the characteristic equation is such that several complex 
roots exist, then the parameter plane curves required to 
realize these roots must all pass through the M-point. If 

I 

the complex roots have the same zeta but different omega 



52 



values, then the constant zeta curve must pass through the 
M-point more than once. In fact, once any point on the 
zeta contour is defined, omega, alpha, and beta are also 
defined, and all roots of the characteristic equation are 
thus fixed. 

2 . Constant Omega Contours 

Within the program, for a given value of omega, zeta 
is varied between zero and one inclusively while omega is 
held constant . 

As with the constant zeta curves, any point on a 
constant omega contour is the omega for a complex root of the 
characteristic equation. By selecting an operating point, 
zeta is also defined whereby a pair of complex conjugate roots 
is established. Again, once any point is chosen on either 
a constant zeta or a constant omega curve, all roots of the 
characteristic equation are established. 

3 . Constant Sigma Contours 

When real roots are to be evaluated, it is 
convenient to return to the characteristic equation: 

™ k 

2 a.S^ = 0 (2-1) 

k=0 ^ 

where, again, a^^ represents a linear combination of constant, 
alpha, and beta terms. S=-<t (a real number) is then an 
equation of a straight line on the parameter plane. If any 
line of constant sigma value passes through the M-point, 
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then the alpha and beta coordinates of the M-point satisfy 
the characteristic equation for a real root located at -a. 

For program considerations, one enters a positive 
value for sigma (corresponding to a real root at and 

the constant sigma contour (straight line) is developed. 

The coordinates of any point on this curve produce a real 
root at -<r. That this is a useful tool, consider that 
system specifications can be achieved by placing a pair of 
complex conjugate roots of the characteristic equation at 
a specific location- To ensure dominance of this root 
pair, the real part of the complex roots so placed should 
be smaller in magnitude than that of the remaining system 
roots. Roots placed at a specified sigma value can thus 
ensure at least one real root whose magnitude is greater 
than the real part of the intended complex conjugate pair. 

4 . Constant Zeta-Omega Curves 

For a fixed value of zeta and omega a pair of complex 
conjugate roots is defined in terms of the expression: 



S = 

The real part of these roots is, thus, defined by the 
zeta-omega product. Note that settling time is defined as 

4 

Ts = . If the product is known, so, too, is the 

duration of the transient response. Thus, by specifying a 
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constant value for the zeta-omega product, any point on the 
contour generated by this value will produce a given 
settling time. 

For any of the parameter plane contours, it is desirable 
to ascertain the values of the characteristic roots for every 
few values of alpha and beta. This feature has been 
incorporated within the tabular output facility as described 
below. 

5 . Tabular Output 

For each of the zeta, omega, and zeta-omega parameter 
plane curves, an arbitrary though reasonable 300 points are 
calculated with which to plot the contours. For the constant 
sigma curves, only two points are needed to define the 
required straight lines (in practice, 4 points are generated 
to ensure that the sigma contours can be plotted within the 
user-defined axes limits). Because of the bulk of data 
points so generated, tabular output is offered as an option 
(as is graphical output), and all points so generated are 
listed for the user. In addition, it is worthwhile to 
calculate system roots for given values of alpha and beta. 
However, computation of roots for each alpha and beta pair 
would cost unnecessary computer time and will likely tax the 
user's patience with the bulk of output so generated. Thus 
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the characteristic roots are generated for every tenth^ 
pair of alpha and beta values. 

6 . Plot Rescaling 

Regardless of the plot scale selected by the user, 
the program generates the full 300 data points for each 
curve requested (4 points for constant sigma lines). When 
the graphical output option is requested, the first family 
of curves is automatically scaled to encompass each and 
every data point. The disadvantage of this technique is 
that , because most activity for the vast majority of systems 
occurs near the physical origin (i.e., alpha=bet a=0 ) , the 
curves may at first appear within only a very small sector 
of the entire plot area, and often they are indistinguishable 
from one another. The advantages of automatic scaling for 
the first set of parameter plane curves far outweigh this 
disadvantage. First, by plotting all available data points, 
the possible limits for alpha and beta are exposed--this is 
important if very large values of alpha and/or beta are 
required to meet the design specifications. Second, for 
some systems the area of activity may not occur near the 
origin, and automatic scaling spares the user the task of 
selecting a sector and possibly missing a sector of interest. 



Although a seemingly arbitrary choice, the generation 
of characteristic roots for every tenth alpha, beta pair 
produces a very tidy output on the commo.nly-used IBM-3278 
computer terminal. 
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Once the user is able to view the panoramic alpha, 
beta parameter curves, it becomes obvious which sector(s) 
are of interest. The user then has the option to rescale 
the set of curves by selecting upper and lower limits for 
the alpha and beta axes. He may continue to rescale the 
family of parameter curves as often as is desirable, and at 
any time the autoscaling option may be recalled. 

The curves generated by the above program are 
sufficient to explore most control system engineering 
problems. The use of these parameter plane contours, and 
their interaction with one another, will be evidenced in 
the next chapter. The source code listing of the program 
is included as Appendix B. 



B. INSTRUCTION TO THE USER 

The parameter plane program is highly interact ive--the 
user is prompted for each required input. A brief 
description of all but the most trivial input items follows 

- Starting value of w ; For most control systems the 
initial value of “ ^is chosen to be zero. Because 

is used in the denominator of certain of the parameter 
plane equations, must be greater than zero. However 
the user may choose arbitrarily close to zero if 
desired . 

- Number of decades: For the majority of control system 

problems a suitable number of decades to be considered 
might be two or three. For higher order systems, it 
would be advisable to start with a slightly larger 
number of decades, especially if the initial “ value 
is small. For subsequent families of curves, ?he 
number of decades can easily be changed. 
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- Constant, alpha, and beta coefficient values: 
Characteristic coefficients are requested from the 
highest to lowest order term. By way of an example, 
a third order characteristic equation might be: 

+ (3a + 3 + 10)S^ + as + (3+5) = 0 

Here, the constant coefficients would be entered in the 
following sequence: 1,10,0,5 while alpha and beta 

coefficients would be entered as 0,3, 1,0 and 0,1, 0,1 
respectively . 

- Zeta values: By convention, values are restricted to 

between zero and one, inclusively. 

- Sigma values: Positive values of sigma correspond to 

negative real roots. Since few, if any, practical 
engineering applications exist for designing a positive 
real root into a system, negative values for sigma are 
disallowed . 

- Omega values: Values for constant omega curves are 

restricted in the lower limit to the starting value, 
and in the upper bound by co^xl0“®^^“®® . 

- Zeta-omega values: As with the constant sigma curves, 

values for constant zeta-omega contours must be greater 
than or equal to zero. 

The user then has the following options: 



1. 


Review 


entries . 


2. 


Change 


any entry. 


3. 


Tabular 


output . 


4. 


Graphical output. 



Remember that tabular output includes 300 data points 
for all but the sigma contours. Characteristic roots are 
displayed for every tenth alpha, beta pair. Because of the 
bulk of output for this option, use it only when necessary. 
If a printed copy of the tabular output is desired, type 
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"record on" before invoking the program. Upon exiting the 
program, type "record off" , after which the user may save 
the preceding terminal session in a listing file designated 
by a name of his choosing. Simply print the listing file, 
which will include all output which has transpired on the 
terminal between the two calls to "record". 

When graphical output is requested, all curves are 
superimposed on the same plot. The first set of curves is 
produced with an autoscaling feature, which plots all 
points calculated (the range of points depends on your 
choice of initial , number of decades, zeta values, omega 
values, etc.). For most characteristic equations only the 
first quadrant (i.e., positive alpha and beta values) will 
be of interest, since negative values usually (but not 
necessarily) imply negative characteristic coefficients and, 
thus, positive roots leading to system instability. The 
nature of the first (autoscaled) set of curves will reveal 
the actual areas of interest for subsequent plots for the 
same system. 

Finally, after each family of curves is plotted, the 
user has six additional options: 



1. 


New 


problem. 


2. 


Same 


problem. 


3. 


Root 


finder . 


4. 


Save 


pro-blem. 
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5. 



Create "DISSPIjA metafile". 



6. Return to main menu. 

Items 1, 4, and 6 are self-explanatory. The remaining 
options deserve some additional explanation. 

- Option 2: This option can be used to re-enter the 

problem at a point prior to actual plotting. Then, 
specific input values can be added or modified, 
tabular and/or graphical output can be requested, and 
entries can be reviewed. It is within this option that 
the graph coordinate axes can be rescaled. If user- 
defined scaling is desired, the minimum and maximum 
values for the axes are requested. 

- Option 3; Although within the tabular output feature a 
set of characteristic roots is produced for every tenth 
pair of alpha and beta values, the bulk of output using 
that option may prove excessive for some applications. 
Here, the user has the option of choosing specific 
values of alpha and beta (e.g. extracted from the 
family of parameter plane curves) and obtaining the 
system roots. 

- Option 5; The program provides a choice from among 
four graphic output devices. Usually the user will 
nominate the TEK618 graphics terminal due to its 
relatively high quality plot resolution. Once the 
user has produced a parameter plane plot to his liking, 
he may wish a final plot of very high resolution. By 
selecting this option, the plot is stored as a DISSPLA 
metafile, and the program is terminated (termination 

of the program is necessary at this point due to an 
anomaly of the DISSPLA graphics package). Simply type 
"DISSPOP" and follow the simple instructions, choosing 
the default options as they are presented. Within the 
"DISSPOP" routine, any of several output devices can be 
called, including the high resolution Versatec plotter 
and the 3800 laser printer. 
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V. PARAMETER PLANE CURVES-GRAPHICAL METHOD 



A. GRAPHICAL SOLUTION 

The algebraic solutions discussed in Chapter III have the 
disadvantage that a fixed value of zeta and omega must first 
be chosen to compute the alpha and beta terms. In some 
instances it is possible to modify the remainder polynomial 
so as to ensure that the specified roots are dominant. 

However, it is not always possible to guarantee that roots 
placed at a specified location can be made dominant. Thus, 
an exhaustive trial-and-error procedure may be required to 
achieve the best values for the various parameters. Trial- 
and-error may also be. required in the design of cascade 
compensators where a specific root location may require 
parameter values that are not physically realizable. In 
these cases, the calculation must be repeated in terms of 
slightly modified specifications; possibly a different means 
of compensation must be used. 

To avoid this trial-and-error analytical approach , one 
can employ the graphical solution. Once a family of curves 
is generated by the program one can, by choosing an M-point 
in the parameter plane, obtain from the curves the n roots 
of the nth order characteristic equation. The trial-and-error 
procedure can then be done visually to reveal an operating 
point which best meets the given specifications. 
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Example 5 . 1 (An Attitude Control System for Large Launch 
Vehicles) 




Figure 5-1 



Figure 5-1 shows the mechanization for a control system for a 
large launch vehicle; Figure 5-2 shows the equivalent block 
diagram 




Figure 5-2 



where the equivalent G(s) is: 



(EqS^ + K^S .+ K^KS^-C) 

G ( s ) “ o p 

+ ^^e^ + S + a C) 

e e 



From G(s) one obtains the system's characteristic equation: 



\ + (l+aQ)S^ + (a C+K )S^ + (K -a C)S^ 

“e ® 

- CK^S - CK^ = 0 

where aQ , a^ , = control system gains 

? ^ = damping ratio of control servo 

^ = natural frequency of control servo 



Gains and are chosen as the system parameters a and 3, 
respectively, to be portrayed in the parameter plane. One 
must then find suitable values for these parameters to yield 
a desired stability margin and a satisfactory transient 
response. For a typical choice of system parameters (such 
as those describing Saturn V), it is assumed that: 

C = 0.717 
e 

“ = 4.71 Hz = 29.594 radians 

e 
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1.0 



C = 0.7 

Then the system's characteristic equation becomes: 

O.OOllS® + 0. 04853^ + 0.53^ + (0.7+a)S^ + (0.35 + B)3^ 

- 0.7“3 - 0.7B = 0 

Various values of aQ and a^ were used to deduce their effect 
on the stability regions, which is indicated in Figure 5-3. 

The numbers of stable and unstable roots, respectively, are 
portrayed in parentheses for each region of the parameter 
plane. 

The analysis precedes as follows: the C =0 contour repre- 

sents boundaries of stability (or relati"ve stability when 
C>0) associated with pairs of complex conjugate roots. The 
3=0 (sigma=0) curve represents real root stability boundaries. 
The region of stability is thus that area bounded by these 
two contours. 3ee Figures 5-4a and 5-4b for a magnified 
view of this area. Any negative alpha-beta pair from within 
this region will exhibit six stable (i.e., all within left- 
half 3-plane) roots and system stability will be assured. 

Note that from the form of the characteristic equation, 
both alpha and beta must be negative to obtain a stable 
system. To illustrate, let us select an arbitrary operating 
point constrained to lie within the lower loop of Figure 5-4b 
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ATTITUDE CONTROL SYSTEM 




ALPHA 



Figure 5-3 

Parameter Plane Curves for 

0.0011 S® + 0.04853^ + 0.53^ + (0.7+a)S^ + (0.3543)5^ 

- 0.73 S - 0.7a = 0 
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BETA 



ATTITUDE CONTROL SYSTEM 




ALPHA 

Figure 5-4a 

p. Parameter Plane Curves for „ 

O.OOllS + 0.04858^ + 0.58'^ + (0.7+a)s3 + (0.35+3)8^ 

- 0.738 - 0.7a = 0 
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BETA 



ATTITUDE CONTROL SYSTEM 




Figure 5-4b 

Parameter Plane Curves for o o 

O.OOllS^ + 0.0485 s 5 + 0.5S4 + (0.7+a)S + (0.35+§)S 

- 0.7BS - 0.7a = 0 
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and also satisfying the requirement that both alpha and beta 
be negative. Superimposed on Figure 5-4b is the constant 
?=0.5 contour, upon which our M-point might be chosen. If 
we select, say, u=-0.15 and B=-0 .02 ( corresponding to 5=0.5 
and cd= 0.5), the system roots are: 

-0.26 + jO.45 

-dominant roots 

-0.26 - jO.45 
-0.32 + jO.12 
-0.32 - j0.12 
-14.6 + jO 
-26.7 + jO 

Coincidentally, the roots associated with our choice of zeta 
and omega are seen to be dominant . The actual choice of an 
operating point may depend on other criteria not discussed 
here . 

Since and (a and 3 , respectively) are functions 
of^i>,5, a^, a^, C,o)^, and 5 ^ » they may be determined for 
various instances of flight by plotting several constant 
zeta and constant omega curves. Actually, and K.^ vary so 
little within the range of values used for C, for specified 
values of ? and SJ, that it becomes possible to choose constant 
values for and . 

This has been a relatively simplistic treatment of a 
complicated control system problem, but it demonstrates the 
power of the parameter plane graphical technique. Knowing 
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nothing more than the system's characteristic equation, the 
static and dynamic stability boundaries may be obtained to 
define the area of overall stability. Of course, other 
system constraints may exist to further limit this area. 

As a note, the root-finding option available within the 
program was used to confirm the numbers of stable and 
unstable roots lying within each region of Figures 5-3 and 5-4. 

Example 5-2 (Alternator Voltage Regulator) 

In designing a voltage regulator of the type shown in Figure 
5-5, we must find values for , and K., that provide 

stability and good transient performance. 



EXCITER ALTERNATOR 




The characteristic equation of the system, including alternator, 
tie-line, etc. is (see reference 8 for details): 



0.00953^ + 0.13253^ + 1.723^ + (K +7.55)3^ + (K +9.1)3 
+ Kq-2.5 = 0 



I 
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We will consider K =0 as fixed and a and B as variable para- 

o 

meters corresponding to K2+7.55 and K^+9.1, respectively. 
Curves of ot versus B are plotted in Figure 5-6 with as 
the variable parameter (for this system, and are 
known to be functions of alone). Since the program plots 
constant zeta curves as a function of varying omega, these 
curves were selected as logical candidates to study the 
problem. 

The beta axis corresponds to the zero value of damping 
constant (sigma) since below the zeta=1.0 curve, the real 
roots (Thaler and Brown 1960) are the negative slopes of the 
tangents drawn from the point in question to the ^=1 curve. 
The machine is then stable for any (a, B) pair between the 
5=0 and S=0 curves. The greatest stability of the machine 
is then possible when both zeta and sigma are largest. 
Further, the best stability can be expected in the region 
bounded by the C=0*3 and 5=1.0 contours (Kabriel 1967). 
Similar stability limits of a and B can be investigated by 
taking as 10, 20, 30 and so on. 

Consider now K^=0 as fixed. Then a and B will represent 
the variable parameters and Kq- 2.5, respectively. 

The characteristic equation then becomes: 

0. 00953^ + 0.1325S‘^ + 1.72S^ + aS^ + 9. IS + B = 0 
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Again, and are known to be functions of alone. Here, 
however, can be solved explicitly [Ref. 8] and one obtains 
three straight-line equations: 

3 = 0 

5.421a -3 = 3.894 
175.63a -3 = 4085.0 

or in terms of system parameters: 

Kq = 2.5 

5.421(K„+7.55)-K = 3.894 - 2.5 

2 o 

175.63(K2+7.55)-Kq = 4085.0 - 2.5 

These are plotted in Figure 5-7. The triangular region 
bounded by these three lines represents a stable region. The 
values of (a , 3 ) within the triangle and hence corresponding 

values of and can be predicted for stable operation. 

The procedure can be repeated for further investigation by 
taking K^=15, 30, 45, etc. 

The analytical parameter plane technique can be used to 
determine the stability limits of Kq and by choosing several 

constant values of . But since has to be selected 

arbitrarily for this purpose, this method becomes cumbersome 
and time-consuming. The method presented here, on the other 
hand, can be used to determine the stable range of either 
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ALTERNATOR VOLTAGE REGULATOR 




Figure 5-6 

Parameter Plane Curves for 
0.00953^ + 0.1325S4 + 1.72S3 + aS^ + BS - 2.5 = 0 
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ALTERNATOR VOLTAGE REGULATOR 




ALPHA 



Figure 5-7 

t- Parameter Plane Curves for 
0.0095S + 0.1325S4 + 1.72S3 + aS^ + 9. IS + S = 0 
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Kq and or and fixing the third parameter. Over and 
above the ability to predict stable operation, the method 
povides a direct measure of the damping at and around a 
chosen operating point. 

Example 5-3 (Lead Compensator) 

Consider the system: 




The characteristic equation is: 

+ (5+P)S^ + (4 + 5P)S^ + (4P+K)S + KZ = 0 

We choose to cancel the pole at S=-l with the zero; thus 
Z=1.0 and the characteristic equation becomes: 

+ (5 + P)S^ + (4 + 5P)S^ + (4P+K)S + K = 0 

Let P =a and K = 6. Then the parameter plane curves are as 
sliown in Figure 5.8, For the coefficients of the character- 
istic equation to remain positive (and thus ensure stability), 
it is convenient to consider only positive values for alpha 
and beta. By inspection we choose ^ = 0.5 and cj^ = 2.0 as a 
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BETA 



LEAD COMPENSATOR 




Figure 5-8 

Parameter Plane Curves for 
S‘^ + (5+ot)S^ + (4+5a)S^ + (4a+6)S + 6 



0 



75 



LEAD COMPENSATOR 




Figure 5-9 

Parameter Plane Curves for 
+ (5+a)S^ + (4+5a)S^ + (4a+3 )s + 3 = 0 
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good operating point, for which a = 4.0 and 3 =24.0. The 
corresponding vicinity of Figure 5-8 is re-scaled in Figure 
5-9. Again, using the root-finding option, the roots 
associated with this (? , pair are shown to be dominant. 

Example 5-4 (Lag Compensator) 

If we are especially concerned with steady-state accuracy 
for a ramp input , it may be advisable to design a lag 
compensator. The parameter plane permits us to consider 
steady-state accuracy while designing the transient response. 
If our system is the same as that considered for the lead 
compensator, the characteristic equation remains: 

+ (5+P)S^ + (4+5P)S^ + (4P+K)S +KZ = 0 

KZ 

But now the error coefficient is K = -r=r . Having three 

e 4P ® 

unknown parameters, K, Z, and P, two must be selected (or 
some combination of two) to be ot and 3 while a numerical value 
is assigned to the third. Once this choice has been made 
the parameter plane curves can be calculated and plotted, and 
the loci of constant can be superimposed. Let Z=0.1, 

P= 0^, and K= 3. The characteristic equation becomes: 

+ (5+oi)S^ + (4+5a)S^ + (4a+3)S + 0.13= 0 

Parameter plane curves are shown in Figures 5-10 and 5-11. 

0 1 ^ 

Lines of ^ ^ = 0.1, 0.2, 0.3, etc. may be superimposed. 

If we select = 0.15 and (similarly to the lead compensator 
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BETA 



LAG COMPENSATOR 




Figure 5-10 

Parameter Plane Curves for 
(5+a)S^ + (4+5a)s2 + (4a+6)S + O.lB = 0 
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LAG COMPENSATOR 




ALPHA 



Figure 5-11 

Parameter Plane Curves for 
+ (5+a)s3 + (4+5a)s2 + (4a+B)S + 0.13 = 0 
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example) ?=0.5 and oj ^=2 . 0 , then a =2 . 14 while 3=12.89. This 

produces dominant roots at -1.0 + jl.7. 

Let us reconsider the systems for which a lead compensator 

and a lag compensator were designed. The characteristic 

equation as well as the error coefficient each contain three 

unknowns (if we assume a numerical value for K ). We can 

e 

imbed the error coefficient in the characteristic equation 

by direct substitution, thereby eliminating one of the 

unknowns. Let us eliminate the gain parameter K - note that 
4PK 

K = 2 • 

Returning to the characteristic equation: 

4 ci o 4PK 

S + (5+P)S + (4+5P)S + (4P+ -„ ■ )S + 4PK = 0 

^ e 

p 

Letting P = a and ^ ~ ^ characteristic equation becomes: 

+ (5+a )S^ + (4+5a )S^ + (4a+4K 3 )S + 4K a =0 

e e 

If a constant value is now chosen for K , say K = 2.0, the 

e e 

characteristic equation becomes: 

+ (5+ct)S^ + (4+5ct)S^ + (4ct+83)S + 8a =0 

for which the parameter plane contours are first shown in 
Figure 5-12 and are further magnified in Figure 5-13. Now 
when any operating point is chosen on the parameter plane 
curve(s), the selected (a, 3) pair generates K^=2.0. Of 
course, this procedure can be repeated for any choice of K^. 
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BETA 



EMBED ERROR COEFFICIENT (KE=2) 




Figure 5-12 

4 Parameter Plane Curves for 

S + (5+0‘)s3 + (4+5a)s2 + (4om-86)S + 8a = o 
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EMBED ERROR COEFFICIENT (KE=2) 




ALPHA 

Figure 5-13 

. Parameter Plane Curves for 

S + (5+a)s3 + (4+5a)s2 + (4a+8S)S + 8a = 0 
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On the parameter plane curves let us choose ==0.5 and 

=2.0 for which: 
n. 

ct = 6.00 = P ■ 

3 = 4.50 = I 
a 

^ = 1.33 = Z 
K = 83 = 36 

KZ 

As a check, = 2.0. The roots associated with the 

selected zeta and omega values are shown to be dominant 
using the root-finding algorithm and are: 

-1.00 + jl.73 

- Dominant Pair 

-1.00 - jl.72 
-1.63 + j 0 
-7.37 + jO 



Example 5-5 

As a final engineering example, consider the system below 
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For this system, the following specifications are to be met: 

1. Set K to the stability limit. 

2. Place a dominant root pair within the following region: 
0.4^C10.7, and 2^co^6. 

3. Both tachometer and acceleration feedback may be used, 
but if possible choose only one. 

From the figure the uncompensated system's open loop transfer 

function is: 

GH = -1 = ^ 

S(S+10)(S +5S+100) 

From which, when expanded, the characteristic equation 
becomes: 

+ 15S^ + 150S^ + lOOOS + K = 0 

To determine the value of K at the stability limit the Fouth 
array is employed: 



1 


150 


K 


15 


1000 


0 


1250 


15K 


0 


1.25xlO^-225K 


0 


0 


15K 


0 


0 
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Here, the stability limit is seen to be K=5555.5. If both 
tachometer and acceleration feedback are used the compensated 
system’s characteristic equation becomes: 

+ 15S^ + (150+5555. 50t)S^ + (1000+5555.53)3 + 5555.5=0 
where a=K , B=K , and K has been set to the stability limit. 

Si L 

Parameter plane curves for this system are shown in Figure 
5-14. 

From these curves, the following analysis can be made. 

The origin of the parameter plane corresponds to the roots 
of the uncompensated system. Since the ?=0 curve passes 
through the origin, two roots are located on the jo) axis of 
the S-plane as was to be expected from the Routh array. The 
remaining two roots are also complex an”d correspond to S=0.8 
and (0=5.0. It is important to note that when an operating 
point involves two different pairs of complex roots, then the 
curves for two different values of omega and two different 
values of zeta must pass through the point. To determine 
which value of omega corresponds to which value of zeta, it 
becomes necessary to refer to the program's tabular data 
output, which is not included here due to lack of space. 

With K =0, the effect of tachometer feedback alone 
corresponds to movement of the M point along the 3 axis. 

In Figure 5-14 the unstable region is determined by an 
inspection of the way the constant zeta curves tend 
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BETA 



TACHOMETER, ACCELERATION FDBK 




Figure 5-14 

Parameter Plane Curves for 

S“^ + 15S^ + (150+5555. 5a )S^ + ( 1000+5555 , 5 B)S + 5555.5 = 0 
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as zeta increases. Since the 3 axis is always in the 
unstable region, it is concluded that tachometer feedback 
alone cannot stabilize the system. 

The effect of acceleration feedback alone can be observed 

by traveling along the a-axis of Figure 5-14. If K is 
varied between 0.01 and 0.06, the system will exhibit two 
pairs of complex roots with the following ranges of values 
for zeta: 

0.3<?<0.5 and 0.25<?<0.32 

If both tachometer and acceleration feedback are used, 
it is seen that tachometer feedback will in general cause 
the zeta of one pair of roots to increase while the zeta 
of the other pair decreases. Acceleration feedback alone 
would appear to be the better choice. 

From the set of curves it is determined that with 
and K^=0.012, four complex roots are located with associated 
values at C=0.45, w=4.0, and ?=0.32, co=13.0. Since 
( 0. 45 ) (4 )=18< < ( 0. 32 ) ( 13 )=4 . 15 , it is apparent that the roots 
at 5=0.45 and |^=4.0 are dominant. The specifications have 
been met and the problem is solved. 

B. MISCELLANEOUS ASPECTS OF THE PARAMETER PLANE 

It can be demonstrated that constant zeta parameter 

plane curves of order two through five originate at a point 
M N 

where alpha = ^ and beta = where M, N, and K are 
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determined by the zero and first power coefficients only. 
Intuition can be used to conclude that constant zeta curves 
of any order originate at this common point which is 
determined only by the zero and first power coefficients. 

An exception is when K=0. In this case the origin point 
depends on higher order coefficients and its location will 
be obvious given a specific problem. If K is not zero, 
the origin is independent of the order of the characteristic 
equation . 

Inspection of the expressions for alpha and beta indicates 
that the shape of the constant zeta curves as omega becomes 
larger is primarily determined by the coefficients of higher 
power, and in general the curves become more complex and 
less well behaved as the order of the characteristic 
equation increases. For a given characteristic equation, 
an increase in complexity can be observed as alpha and beta 
appear in more coefficients. 

All constant zeta curves tend to plus or minus infinity. 
The relative magnitudes of the coefficients determine whether 
the limit is plus or minus infinity. It is therefore 
necessary to choose a frequency range of interest before 
plotting the curves, thus limiting the analysis to one 
"window" of the infinite plane. 

^ Since no stability criteria, either relative or 
absolute, has been established for the parameter plane, it 
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is necessary to base the stability analysis on observing 
which way the curves tend as omega and zeta are varied. 

For this reason it is worthwhile to plot curves for as many 
values of zeta, omega, sigma, and if desired, zeta-omega, 
as are necessary to ascertain the pattern. 
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VI. CONCLUSIONS 



Parameter plane techniques have been applied to the 
compensation of linear control systems. General equations 
have been derived for the cases of feedback, cascade, and 
combination feedback-cascade compensation, to enable one to 
place a pair of complex conjugate roots at a specific 
location in the S-plane, while simultaneously satisfying 
the steady-state accuracy requirements. A dominancy 
technique has been introduced whereby once a pair of complex 
roots is fixed, the remaining roots of the characteristic 
equation can be manipulated to ensure that the specified 
roots are dominant. 

The impetus for development of a parameter plane program 
was to provide the user with a quick, simple means of 
obtaining the information available in the analytical 
solution to control system compensation, while avoiding the 
painstaking labor of trial-and-error analysis inherent in 
that technique. Several, practical engineering examples 
have been presented to demonstrate the superiority of the 
graphical technique. To date, no other package is known to 
offer the fully interactive and comprehensive capabilities 
of the parameter plane program. 
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By itself the program allows one to design a control 
system compensation model for most systems. However, for 
some lightly damped systems containing mechanical resonances, 
the amount by which zeta or omega are incremented in the 
parameter plane equations may be so large as to not detect 
the resonance peaks. This information would be available from 
either a root-locus or Bode analysis. For still other 
systems, one might be interested in the way the roots of the 
characteristic equation extend from the open loop poles and/or 
zeros. Since the parameter plane equations are calculated 
using only the characteristic equation, no knowledge of open 
loop poles or zeros is available. Again, a root-locus method 
would reveal this information. Incorporated into one 
comprehensive package which includes Bode and root-locus 
analyses, the program provides the capability to investigate 
the entire gamut of linear control system architecture. 

A basis for further investigation involves plotting the 
parameter plane contours for systems that are non-linear 
in the alpha and beta terms — i.e., those systems which contain 
alpha-beta product terms. Although the recursion technique 
used in this text has distinct advantages over the matrix 
approach for the linear case, as pointed out earlier, the 
matrix technique would be the method of choice for the non- 
linear case. 
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APPENDIX A 
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APPENDIX B 



PARAMETER PLANE PROGRAM 



SUBROUTINE LPARAM 
DIMENSION AC350), B(350), 

X AGC350), BGC350), AJ(IOO), BJ(IOO), CJ(IOO), 

X ZETA(IOO), SIGMA(IOO), W(lOO), ZW(IOO) 

C 

C INITIAL ASSIGNMENTS 

CHARACTERXA SCHAR/'S V,YES/'Y ’/,NOO/’N BLANK/’ ’/, 

X ENCHAR/'E •/, WNCHAR/ 'WN '/, NDCHAR/ ’ ND •/, NOCHAR/ ’ NO '/ , 

X A JCHAR/ ’ AJ ’/ , BJCHAR/ ' B J ’/ , CJCHAR/ ' CJ '/ , NSCHAR/ ’ NS '/ , 

X NZCHAR/’NZ’/,ZHCHAR/'ZW’/,NWCHAR/'NW'/,PRCHAR/’PR'/, 

X NCCHAR/'NC’/,TICHAR/'TI'/ 

CHARACTERXA TABLE, GRAPH, CHANGE, REPLY, OPT, LABEL(9) 

COMMON /SAVE/ LABEL, WN, ND, N02, NC, CJ, AJ, BJ, 

X NZ, ZETA, NS, SIGMA, NW, W, NZW, ZW, 

X XMIN, XMAX, YMIN, YMAX 

C 

C DATA ENTRY FROM FILE OR CONSOLE? 

100 MINMAX = 1 
GRD = 0. 

CHANGE = BLANK 

CALL EXCMS( ’CLRSCRN’) 

WRITE(6,500) 

CALL READC (REPLY) 

IF (REPLY .NE. ’D') GOTO lOI 
CALL GETIT 
MINMAX = 0 
GO TO 200 

101 CONTINUE 
C 

C GET CURVE TITLE 

102 CONTINUE 

CALL EXCMS( ’CLRSCRN’) 

WRITE(6,50I) 

CALL READL (LABEL) 

CALL ASTER ( LABEL , LABEL ) 

IF ( CHANGE .EQ. TICHAR ) GO TO 200 
C 

C GET STARTING VALUE OF WN 

103 CONTINUE 
WRITE(6,502) 

CALL READR (WN) 

IF (WN) lOA, 104,105 

104 WRITE(6,503) 

GO TO 103 

105 IF ( CHANGE .EQ. WNCHAR) GOTO 200 
C 

C GET THE NUMBER OF DECADES CONSIDERED 

106 CONTINUE 
WRITE(6,504) 

CALL READI (ND) 

IF ( CHANGE .EQ. NDCHAR ) GOTO 200 
C 

C GET THE ORDER OF THE CHARACTERISTIC EQN 

107 CONTINUE 
WRITE(6,505) 

CALL READI (N02) 

NC = N02+1 

IF (CHANGE .EQ. NOCHAR ) GOTO 200 
C 
C 
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noon on on no on on on on on 



GET THE NUMBER OF CONSTANT ZETA CURVES 

108 CONTINUE 

CALL EXCMS C’CLRSCRN') 

NRITE(6,506) 

CALL READI (NZ) 

IF (NZ .LT. 1) GOTO 110 

GET THE VALUES OF ZETA 

NRITE(6,507) 

DO 110 I = 1,NZ 

109 WRITEC6, 508)1 

CALL READR (ZETA(D) 

IF (ZETA(I) .LT. 0. .OR, ZETA(I) .GT. 1.) NRITE(6,509) 

IF (ZETA(I) .LT. 0. .OR. ZETA(I) ,GT, 1.) GO TO 109 

110 CONTINUE 

IF ( CHANGE. EQ. NZCHAR ) GOTO 200 

GET THE NUMBER OF CONSTANT SIGMA CURVES 

111 CONTINUE 

CALL EXCMS (»CLRSCRN*) 

NRITE(6,510) 

CALL READI (NS) 

IF (NS .LT. 1) GOTO 113 

GET THE VALUES OF SIGMA 

DO 113 I = 1,NS 

112 NRITE(6,511) I 

CALL READR (SIGMA(D) 

IF (SIGMA(I) .LT. 0.) WRITE (6,512) 

IF (SIGMA(I) .LT. 0.) GO TO 112 

113 CONTINUE 

IF ( CHANGE .EQ. NSCHAR ) GOTO 200 

GET THE NUMBER OF CONSTANT WN CURVES 

llA CONTINUE 

CALL EXCMS (»CLRSCRN*) 

WRITE(6,513) 

CALL READI (NW) 

IF (NW .LT. 1) GOTO 116 

GET THE WN VALUES 

WNMAX = WNXIOXXND 
DO 116 I = 1,NW 

115 WRITE(6,514) I 
CALL READR (W(D) 

IF (W(I) .LT. WN .OR, W(I) ,GT. WNMAX) WRITE (6,515) WN, WNMAX 
IF (W(I) .LT. WN .OR. W(I) ,GT. WNMAX) GO TO 115 

116 CONTINUE 

IF ( CHANGE ,EQ. NWCHAR ) GOTO 200 

GET THE NUMBER OF CONSTANT ZETAxWN CURVES 

117 CONTINUE 

CALL EXCMS CCLRSCRN*) 

WRITE(6,516) 

CALL READI (NZW) 

IF (NZW .LT. 1) GOTO 119 

GET THE ZXWN VALUES 

DO 119 I = 1,NZW 

118 WRITE(6,517) I 
CALL READR (ZW(D) 

IF (ZW(I) .LE. 0.) WRITE (6,518) 

IF (ZW(I) .LE. 0.) GO TO 118 

119 CONTINUE 

' IF ( CHANGE .EQ. ZWCHAR ) GOTO 200 
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120 

121 

122 

125 

124 

125 



GET CONSTANT COEFFICIENTS 

CONTINUE 

CALL EXCMS CCLRSCRN’) 

WRITE(6,519) 

DO 121 I = NC,1,-1 
II = I-l 

WRITE(6,520) II, I 
CALL READR (CJ(D) 

CONTINUE 

IF ( CHANGE .EQ. CJCHAR ) GOTO 200 

GET ALPHA COEFFICIENTS 

CONTINUE 

CALL EXCMS C’CLRSCRN’) 

WRITE(6,521) 

DO 123 I = NC,1,-1 
II = I-l 

WRITE(6,522) II, I 
CALL READR (AJ(D) 

CONTINUE 

IF (CHANGE .EQ. AJCHAR ) GOTO 200 

GET BETA COEFFICIENTS 

CONTINUE 

CALL EXCMS (’CLRSCRN’) 

WRITE(6,523) 

DO 125 I = NC,1,-1 
II = I-l 

WRITE(6,524) II, I 
CALL READR (BJ(D) 

CONTINUE 
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REVIEN ENTRIES 

200 CONTINUE 

CALL EXCMS CCLRSCRN') 

WRITE(6,525) 

CALL READC (REPLY) 

IF (REPLY .NE. YES) GOTO 209 
CALL EXCMS ('CLRSCRN') 

WRITE(6,526) 

WRITE(6,527) (LABEL(I), 1=1, 9) 
WRITE(6,528) 

WRITE(6,529) ND, N02, NZ, NS, NW, NZW 

WRITE(6,530) WN 

WRITE(6,531) 

IF (NZ .LT. 1) GOTO 201 

NRITE(6,532) (ZETA(M) , M=1 , NZ) 

GOTO 202 

201 WRITE(6,533) 

202 CONTINUE 
WRITE(6,534) 

IF (NS .LT. 1) GOTO 203 

NRITE(6,532) ( SIGMA(M) , M=1 , NS) 

GOTO 204 

203 NRITE(6,533) 

204 CONTINUE 
NRITE(6,535) 

IF (NW .LT. 1) GOTO 205 

WRITE(6,532) ( W(M) ,M=1 , NW) 

GOTO 206 
WRITE(6,533) 

206 CONTINUE 

WRITE(6,536) 

IF (NZW .LT. 1) GOTO 207 

WRITE(6,532) (ZW(M) ,M=1 , NZW) 

GOTO 208 
WRITE(6,533) 

208 CONTINUE 
WRITE (6,576) 

CALL EXCMS( »CLRSCRN» ) 

WRITE(6,537) 

WRITE(6,532) (CJ ( N) , N=NC, 1 , -1 ) 
WRITE(6,538) 

WRITE(6,532) ( AJ ( N) , N=NC, 1 , -1 ) 
WRITE(6,539) 

WRITE(6,532) ( BJ ( N) , N=NC, 1 , -1 ) 
WRITE(6,540) 

WRITE(6,541) XMIN, XMAX, YMIN, YMAX 

209 CONTINUE 
WRITE(6,542) 

CALL READC (CHANGE) 

IF (CHANGE .NE. YES) GOTO 210 
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OCJOO CJ OOCJCJOCJCJOOCJOOOOCJUOCJCJCJOCJCJOCJCJCJOOCJCJOCJCJOCJCJCJ 



CHANGE ROUTINE 



CALL EXCMS CCLRSCRN') 
WRITE(6,543) 



CAL 


L READC 


(CHANGE) 






IF 


(CHANGE 


.EQ. 


TICHAR) 


GOTO 


102 


IF 


(CHANGE 


.EQ. 


NNCHAR) 


GOTO 


103 


IF 


(CHANGE 


.EQ. 


NDCHAR) 


GOTO 


106 


IF 


(CHANGE 


.EQ. 


NOCHAR) 


GOTO 


107 


IF 


(CHANGE 


.EQ. 


NZCHAR) 


GOTO 


108 


IF 


(CHANGE 


.EQ. 


NSCHAR) 


GOTO 


111 


IF 


(CHANGE 


.EQ. 


NNCHAR) 


GOTO 


11^ 


IF 


(CHANGE 


.EQ. 


ZNCHAR) 


GOTO 


117 


IF 


(CHANGE 


.EQ. 


CJCHAR) 


GOTO 


120 


IF 


(CHANGE 


.EQ. 


AJCHAR) 


GOTO 


122 


IF 


(CHANGE 


.EQ. 


BJCHAR) 


GOTO 


12^ 


IF 


(CHANGE 


.EQ. 


ENCHAR) 


GOTO 


^05 


IF 


(CHANGE 


.EQ. 


PRCHAR) 


GOTO 


100 


IF 


(CHANGE 


.EQ. 


NCCHAR) 


GOTO 


210 



GOTO 209 

210 CONTINUE 
WRITE(6,5A4) 

CALL READC (TABLE) 
WRITE(6,5A5) 

CALL READC (GRAPH) 

211 IF (NZ) 213,213,212 

212 G = (10.XXND)**(l./299.) 

213 CONTINUE 
2 = l.OE-60 



97 



c 

c 

300 



C 

C 



C 



C 



301 

302 



303 
C 

304 



C 



305 

306 
C 

C 



C 

307 



C 

C 

C 

C 

C 

C 



308 

309 



X 

X 

X 



CONSTANT ZETA PLOTS 

IF (NZ) 309,309,300 

CALL EXCMS (»CLRSCRN») 

IF (TABLE .EQ. YES .OR. GRAPH . EQ . YES) WRITE (6,546) 



DO 308 M=1,NZ 

IF (TABLE .EQ. YES) WRITE (6,547) 

J = 0 
R = 0. 

WNA == WN 

DO 306 L=l,300 
D1 = 0.0 
D2 = 0.0 
Cl = 0.0 
C2 = 0.0 
B1 = 0.0 
B2 = 0.0 

DO 303 N=1,NC 
K = N-1 

IF (K) 302,301,302 
U = 0.0 
U1 = -1.0 

U2 = 2.05^ZETA(M)3^U-U1 
D1 = (-1.0)XXKXCJ(N)5^WNA9^5^K3(U1 + D1 
D2 = (-1.0)5^XKXCJ(N)XWNAXXKXU+D2 
Cl = (-1.0)5(XKXBJ(N)3^WNAX3^K3^U1+C1 
C2 = (-1.0)XXKXBJ(N)5^WNA)(XK3(U+C2 
B1 = (-1.0)XXKXAJ(N)XWNAX3(KXU1 + B1 
B2 = (-1 .0)XXK3(AJ(N)XWNA3(5(KXU+B2 
U1 = U 
U = U2 
CONTINUE 

IF (ABS(B13(C2-B2XC1)-Z) 306,306,304 
J = J + 1 
R = R+1. 

A(J) = (C1XD2-C23^D1)/(B13(C2-B23^C1) 

B(J) = (B25(D1-B13^D2)/(B19^C2-B23(C1) 

IF (TABLE .NE. YES) GO TO 306 

WRITE(6,548) A(J), B(J), WNA, ZETA(M) 

IF (R/10. - J/10) 306,305,306 

CALL ROOTS (A(J), B(J), A J , BJ, CJ, N02) 

CALL EXCMS(»CLRSCRN') 

WRITE (6,547) 

WNA = G5^WNA 

CALL EXCMS(»CLRSCRN») 

IF (J .GT. 0) GOTO 307 
WRITE(6,549) 

GOTO 308 

IF (GRAPH .EQ. YES) CALL PLOTD( A, B , J , . FALSE . , 

LABEL, »ALPHA$», 'BETA$', 
MINMAX, » Z=$»,ZETA(M), 
XMIN,XMAX,YMIN,YMAX,GRD) 

IF (GRAPH .EQ. YES) GRD = GRD+1. 

CONTINUE 

CONTINUE 
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CONSTANT SIGMA PLOTS 

IF (NS) 335,335,325 
325 CALL EXCMS( 'CLRSCRN' ) 

IF (TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,550) 



XINC = (XMAX - XMIN)/299. 

YINC = (YMAX - YMIN)/299. 

DO 33A M=1,NS 
XPT = XMIN 
YPT = YMIN 
DD = CJ(1) 

CC = BJ(1) 

BB = AJ(1) 

J = 0 
R = 0. 



DO 326 N=2,NC 
K = N-1 

DUMMY5=SIGMA(M)XXK 
DD = (-1.0)XXKXCJ(N)XDUMMY5+DD 
CC = (-1.0)XXKXBJ(N)XDUMMY5+CC 
BB = (-1.0)*XKXAJ(N)XDUMMY5+BB 

326 CONTINUE 

IF (CC .EQ. 0. .AND. BB .EQ. 0.) GOTO 334 
IF (CC) 327,327,330 

327 DO 329 L=l,300 

J = J + 1 
R = R+1. 

A(J) = XPT 

B(J) = (-BBXA(J)-DD)/CC 
IF (TABLE .NE. YES) GOTO 329 
WRITE(6,552) A(J), B(J), SIGMA(M) 

IF (R/10. - J/10) 329,328,329 

328 CALL ROOTS (A(J), B(J), AJ, BJ, CJ, N02) 
CALL EXCMS ('CLRSCRN') 

WRITE(6,550) 

329 XPT = XPT + XINC 
GO TO 333 



330 DO 332 L=l,300 

J = J+1 
R = R+1. 

B(J) = YPT 

A(J) = (-CC5(B(J)-DD)/BB 
IF (TABLE .NE. YES) GOTO 332 
WRITE(6,552) A(J), B(J), SIGMA(M) 

IF (R/10. - J/10) 332,331,332 

331 CALL ROOTS (A(J), B(J), AJ, BJ, CJ, N02) 
CALL EXCMS ('CLRSCRN') 

•WRITE(6,550) 

332 YPT = YPT + XINC 



333 



X 

X 

X 



CALL EXCMS( 'CLRSCRN' ) 

IF (GRAPH .EQ. YES) CALL PL0TD( A, B, J , . FALSE . , 

LABEL, 'ALPHA$', 'BETA$', 
MINMAX,' S=$',SIGMA(M), 
XMIN, XMAX, YMIN, YMAX, GRD) 



IF (GRAPH .EQ. YES) GRD = GRD+1 . 



334 CONTINUE 



335 CONTINUE 
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c 

C 

350 

C 

C 

C 



C 



C 



351 

352 



353 

C 

35‘i 



C 



355 



356 

C 

C 



C 

357 

X 

X 

X 

C 

C 

358 
C 

359 
C 



CONSTANT ZETA-OMEGA PLOTS 
IF (NZW) 359,359,350 
CALL EXCMS( ’CLRSCRN’) 

IF (TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,553) 

XWN = WN 

DO 358 M=1,NZW 

IF (TABLE .EQ. YES) WRITE (6,55A) 

J = 0 
R = 0. 

AZETA = 1./300. 

DO 356 L=l,300 

XWN = ZW(M)/AZETA 

D1 = 0.0 

D2 = 0.0 

Cl = 0.0 

C2 = 0.0 

B1 = 0.0 

B2 = 0.0 

DO 353 N=1,NC 
K = N-1 

IF (K) 352,351,352 
Q1 = 0.0 
Q = -1.0/XWNXX2 
D2 = CJ(N)XQ1+D2 
C2 = BJ(N)XQ1+C2 
B2 = AJ(N)xq 1+B2 
D1 = CJ(N)XQ+D1 
Cl = BJ(N)XQ+C1 
B1 = AJ(N)xq+b1 
Q2 = -2.0XZW(M)XQ1-XWNXX2XQ 
Q = Q1 
Q1 = Q2 

IF (ABS(B1XC2-B2XC1)-Z) 356,356,354 
J = J+1 
R = R+1. 

A(J) = (C1XD2-C2XD1)/(B1XC2-B2XC1) 

B(J) = (B2XD1-B1XD2)/(B13(C2-B2XC1) 

IF (TABLE .NE. YES) GO TO 356 
WRITE(6,552) A(J), B(J), ZW(M) 

IF (R/10. - J/10) 356,355,356 

CALL ROOTS (A(J), B(J), AJ , BJ, CJ, N02) 

CALL EXCMS( 'CLRSCRN’) 

WRITE (6,554) 

AZETA = AZETA+(1 ./300. ) 

CALL EXCMS( 'CLRSCRN' ) 

IF (J .GT. 0) GOTO 357 
WRITE(6,549) 

GOTO 358 

IF (GRAPH .EQ. YES) CALL PLOTD( A, B, J , . FALSE . , 

LABEL, 'ALPHAS', 'BETA$', 
MINMAX, 'ZW=$',ZW(M), 
XMIN,XMAX,YMIN,YMAX,GRD) 

IF (GRAPH .EQ. YES) GRD = GRD+1 . 

CONTINUE 

CONTINUE 
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CONSTANT OMEGA PLOTS 

375 IF (NN) 385,385,376 

376 CALL EXCMS('CLRSCRN') 

IF (TABLE .EQ. YES .OR. GRAPH . EQ . YES) WRITE (6,555) 



DO 38A M=1,NW 

IF (TABLE .EQ. YES) WRITE (6,547) 

J = 0 
R = 0. 

AZETA =0.0 

DO 382 L=l,300 
D1 = 0.0 
D2 = 0.0 
Cl = 0.0 
C2 = 0.0 
B1 = 0.0 
B2 = 0.0 

DO 379 N=1,NC 
K = N-1 

IF (K) 378,377,378 

377 U = 0.0 

U1 = -1.0 

378 U2 = 2.03CAZETAXU-U1 

D1 = (-1.0)XXK3^CJ(N)XW(M)9(3^KXU1 + D1 
D2 = (-1.0)9«XK3(CJ(N)XW(M)XXK3(U+D2 
Cl = (-1.0)3«XK5(BJ(N)XW(M)5(XKXU1+C1 
C2 = (-1.0 )XXKxbJ(N)5«W(M)3(3«KxU+C2 
B1 = (-1.0)XXKXAJ(N)XW(M)3«3(KXU1 + B1 
B2 = (-1.0)5«9(K3(AJ(N)xW(M)5(9(KxU+B2 
U1 = U 
U = U2 

379 CONTINUE 

IF (ABS(B19^C2-B2)€C1)-Z) 382,382,380 

380 J = J+1 

R = R+1. 

A(J) = (C1xD2-C2XD1)/(B1XC2-B2XC1) 
B(J) = (B23(D1-B1XD2)/(B15(C2-B23(C1) 

IF (TABLE .NE. YES) GO TO 382 
WRITE(6,548) A(J), B(J), W(M), AZETA 
IF (R/10. - J/10) 382,381,382 

381 CALL ROOTS (A(J), B(J), A J , BJ, CJ, N02) 
CALL EXCMS( 'CLRSCRN*) 

WRITE (6,547) 

382 AZETA = AZETA+( 1 . /299 . ) 

CALL EXCMS( 'CLRSCRN') 

IF (J .GT. 0) GOTO 383 
WRITE(6,549) 

GOTO 384 



383 

X 

X 

X 



IF (GRAPH .EQ. YES) CALL PLOTD( A, B , J , . FALSE . , 

LABEL, 'ALPHAS*, *BETA$*, 
MINMAX, ' W=$»,W(M), 

XMI N , XMAX , YMI N , YMAX , GRD ) 



384 CONTINUE 



385 CONTINUE 
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400 



C 

C 



401 



C 

C 



C 

C 



C 

C 



C 

c 



IF (GRAPH .EQ. YES) CALL PLOTDC 0 . , 0 . , 0 . . . TRUE . , 
X » 

X 0. , » $S-9.7531, 

X 0.,0.,0.,0.,GRD) 

IF (GRAPH .EQ. YES) GRD = GRD+1 . 

GRD = 0. 

CALL EXCMS( »CLRSCRN») 

IF (OPT .EQ. YES) CALL DONEPL 
IF (OPT .EQ. YES) STOP 
WRITE(6,557) 

NRITE(6,556) 

WRITE(6,557) 

WRITE(6,558) 

WRITE(6,559) 

WRITE(6,560) 

WRITE(6,561) 

NRITE(6,562) 

WRITE(6,563) 

WRITE(6,557) 

CALL READI (IANS) 

IF (IANS .GT. 6 .OR. IANS .LT. 1) GOTO 400 
GO TO (100, 404, 401, 402, 403, 405) IANS 

ROOT FINDER OPTION 

CALL EXCMS(*CLRSCRN») 

WRITE(6,564) 

WRITE(6,565) 

CALL READR (ALPHA) 

WRITE(6,566) 

CALL READR (BETA) 

CALL ROOTS (ALPHA, BETA, AJ, BJ, CJ, N02) 

GO TO 400 



402 



CALL SAVIT 
GO TO 400 



SAVE OPTION 



403 



CREATE DISSPLA METAFILE OPTION 

CALL EXCMS( »CLRSCRN») 

WRITE(6,567) 

WRITE(6,568) 

WRITE(6,569) 

CALL READC (OPT) 

IF (OPT .NE. YES) GOTO 400 
CALL EXCMS( »CLRSCRN») 

WRITE(6,570) 

CALL DONEPL 



CALL META 
GO TO 211 



COMPRS = SUBROUTINE TO LET DONEPL FINISH 



404 



405 



SAME PROBLEM OPTION 

WRITE(6,571) 

CALL READI (MINMAX) 

IF (MINMAX .EQ. 1) GO TO 200 
CALL EXCMS( »CLRSCRN*) 

WRITE(6,572) 

CALL READR (XMIN) 

WRITE(6,573) 

CALL READR (XMAX) 

WRITE(6,574) 

CALL READR (YMIN) 

WRITE(6,575) 

CALL READR (YMAX) 

GO TO 200 

CONTINUE 

RETURN 
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500 



501 

502 

503 
50 ^ 

505 

506 

507 

508 

509 

510 

511 

512 

513 
51 ^ 
515 



516 

517 

518 

519 

520 

521 

522 

523 
52 ^ 

525 

526 

527 

528 

529 

530 

531 

532 

533 

534 

535 

536 

537 

538 

539 

540 

541 

542 

543 



FORMATS STATEMENTS 

F0RMATC5X, 'THIS IS THE INTERACTIVE PARAMETER PLANE PROGRAM...',/, 
X 5X, 'THE USER WILL BE PROMPTED FOR VARIOUS INPUTS.',////, 

X 5X,'WILL YOU BE ENTERING DATA FROM A CONSOLE OR DATAFILE?' 

X / 15X '"D” OR *'C"') 

FORMATC///, '"enter TITLE TO APPEAR FOR THIS FAMILY OF CURVES') 
F0RMAT(/,1X, 'WHAT IS THE STARTING VALUE OF OMEGAN (WN>0.0)?') 
F0RMAT(/,1X, 'WN MUST BE GREATER THAN ZERO - TRY AGAIN') 
F0RMAT(/,1X, 'HOW MANY DECADES PAST WN ARE DESIRED? (ND)') 
FORMATC/,' WHAT IS THE ORDER OF THE CHARACTERISTIC EQUATION? (NO ) ' 
F0RMAT(/,1X, 'HOW MANY CONSTANT ZETA CURVES ARE DESIRED? (NZ)') 
FORMATC/,' ENTER THE VALUES OF ZETA TO BE USED IN COMPUTATION...' 
F0RMATC5X, 'ZETAC ',12, ') = ?') 



F0RMATC5X, 
FORMATC/, ' 
X (NS)') 
F0RMATC5X, 
F0RMATC5X, 



'ZETA MUST LIE BETWEEN 0 
HOW MANY CONSTANT SIGMA 



AND 1, INCLUSIVE - TRY AGAIN') 
(REAL ROOT) CURVES ARE DESIRED? 



REAL ROOT - TRY ANOTHER') 
ARE DESIRED? CNW)') 



'SIGMAC ',12, ') = ?') 

'NEGATIVE SIGMA MEANS POSITIVE 
FORMATC/, IX, 'HOW MANY CONSTANT WN CURVES 
F0RMATC5X, 'W(',I2, ') = ?') 

FORMAT C5X,'WN NOT WITHIN PLOTTABLE RANGE', 

X /,5X,'Y0UR USABLE RANGE IS' 

X /,10X,F10.2, ' TO ',F10.2) 

FORMATC/, IX, 'HOW MANY CONSTANT Z^WN CURVES ARE DESIRED? CNZW)') 
F0RMATC5X, 'ZWC ',12, ') = ?') 

F0RMATC5X, 'NON-POSITIVE Z-WN MEANS POSITIVE ROOT - TRY ANOTHER') 
FORMATC/, IX, 'ENTER THE CONSTANT COEFFICIENTS...') 

F0RMATC5X, '— SXXM2, '— CJC',12,') = ?') 

FORMATC/, IX, 'ENTER THE ALPHA COEFFICIENTS...') 

F0RMATC5X, '— SX5^M2, '— AJC',12,') = ?') 

FORMATC/, IX, 'ENTER THE BETA COEFFICIENTS...') 

F0RMATC5X, '— S3e5^',I2, ' BJC',12,') = ?') 

FORMATC/,' WANT TO REVIEW YOUR ENTRIES BEFORE RUNNING? 

FORMAT C/,10X, ' GRAPH TITLE') 

(1X,9A4) 

C/,8X,2HND,8X,2HN0,8X,2HNZ,8X,2HNS,8X,2HNW,7X,3HNZW) 
(6110) 

C/,10X, 'INITIAL VALUE OF OMEGA 
C/,10X, ' ZETA ') 

C8E10.3) 

CIX, ' NO VALUE ' ) 

C/,10X, 'SIGMA ') 

W ') 

'ZW ' ) 

'CONSTANT COEFFICIENTS 
'ALPHA COEFFICIENTS 

'BETA COEFFICIENTS 



CY/N)') 



FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMATC/, 
FORMATC/, 



',F10.5) 



C/,10X, 

C/,10X, 

C/,10X, 

C/,10X, 

C/,10X, 

C/,10X,'XMIN XMAX YMIN YMAX') 
(1X,4E10.3) 

WANT TO MAKE ANY CHANGES? (Y/N)') 
WHAT VARIABLE/AREA DO YOU WISH TO 



IN 

IN 

IN 



DECENDING 

DECENDING 

DECENDING 



ORDER' 

ORDER' 

ORDER' 



X5X, 'TITLE TI 



OMEGA START.. WN 



CHANGE?',//, 
t DECADES ND',/ 
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544 

545 

546 

547 

548 

549 

550 

551 

552 

553 

554 

555 

556 

557 

558 

559 
56 0 

561 

562 

563 

564 

565 

566 

567 



568 



569 



X5X, 

X5X, 

X5X, 

X5X, 

XIX, 



ORDER NO 

WN 

TERMS, 



CONST SIGMA. .NS*,/ 



BETA TERMS. 
NO CHANGE.. 



.BJ* ,/ 
.NC*,///, 



CONST ZETA. . . NZ 

'CONST WN NW CONST Z-WN . . . ZW* ,// 

'CONST TERMS.. CJ ALPHA TERMS.. AJ 

'END E NEW PROBLEM. .PR 

'ENTER TWO DIGIT CODE. . . *) 

FORMATC/,* DO YOU WANT TABULATED DATA ON 
FORMATC/,* DO YOU WANT THE GRAPHS ON THE 
FORMAT (1H1,10X, 'CONSTANT ZETA CURVES') 

FORMAT(/,10X, 'ALPHA BETA OMEGA ZETA') 

FORMAT (4E16.5) 

FORMATC/,' DUE TO PLOT RESTRICTIONS, A COMPLETE GRAPH CANNOT BE OU 
XTPUT. ' ) 

FORMAT (1H1,10X, 'CONSTANT 
(/,10X, 'ALPHA 
(3E16.5) 

(1H1,10X, 'CONSTANT 
(/,10X, 'ALPHA 
(1H1,10X, 'CONSTANT 
OPTION NO. I 



FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
F0RMATC5X, 
F0RMATC5X, 
F0RMATC5X, 
F0RMATC5X, 
FORMATC5X, 
F0RMATC5X, 
F0RMATC5X, 
F0RMATC5X, 
F0RMATC5X, 

Y ^Y 

FORMATC/// 
FORMATC/, 5X, 
FORMATC5X, 

X 5X, 

X 5X, 

X 5X, 

X 5X 

FORMATCSX' 

X 5X, 

X 5X 

FORMATC/// 

X 15X 



SIGMA CURVES') 
BETA 



THE SCREEN? CY/N)') 
TERMINAL? CY/N)') 

OMEGA 



SIGMA') 



ZETA-OMEGA CURVES') 

BETA ZETA-OMEGA' 

OMEGA CURVES') 

OPTION 



'/ 

1 



1 I NEW PROBLEM 

2 I SAME PROBLEM 

3 I ROOT FINDER 

4 I SAVE DATA 

5 I SAVE GRAPH IN DISSPLA METAFILE 

6 I RETURN TO MAIN MENU 

ENTER AN ALPHA-BETA PAIR, AND THE ROOTS OF YOUR 
SYSTEMS CHARACTERISTIC EQUATION WILL BE RETURNED 

E ALPHA VALUE 
THE BETA VALUE 

YOU NOW HAVE THE OPTION OF STORING THE LAST SET OF 
CURVES IN A DISSPLA METAFILE. THIS ALLOWS RETRIEVAL'/, 
OF DATA AT A LATER TIME FOR ROUTING TO ANY OF '/, 

SEVERAL OUTPUT DEVICES CTEK618, 3800 LASER PRINTER, '/, 
VERSATEC PLOTTER, ETC.) 

IF YOU CHOOSE THIS OPTION, THE PROGRAM MUST 
TERMINATED - THIS CANNOT BE AVOIDED WITHOUT 
CATASTROPHIC RESULTS. 

5X,'D0 YOU WISH TO USE THIS OPTION? 

'"Y" OR "N" 



5X, 'ENTER 
ENTER 






BE 



570 FORMATC/////, 5X, 'IF YOU WISH GRAPHIC OUTPUT, TYPE: 



571 

572 
57 3 

574 

575 

576 



X 15X, '"DISSPOP" 

X 5X,'AND FOLLOW THE INSTRUCTIONS . . . 

FORMAT C///,' AUTOSCALE OR USER-DEFINED LIMITS FOR CURVES? 
X/,' l=AUTOSCALE; 0=USER-DEFINED ' ) 



'//) 

'/, 

'/. 

V//) 

'/, 

'/) 

'/. 

'/ 

') 



FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

END 



C//,' INPUT MINIMUM VALUE FOR X CX-MIN)') 
C/,' INPUT MAXIMUM VALUE FOR X CX-MAX)') 

C/,' INPUT MINIMUM VALUE FOR Y CY-MIN)') 

C/,' INPUT MAXIMUM VALUE FOR Y CY-MAX) ' ) 

) 
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SUBROUTINE PLOTD — GRAPHS WILL BE PRODUCED ON UPRIGHT II X lA 
PAGE WITH 9 INCH AXES. IF USER SELECTS 'AUTOSCALE' FEATURE, 
SUBROUTINE PLDOlO (INTERNAL TO PLOTD) FINDS MIN AND MAX FOR EACH = 
AXIS AND SCALES ACCORDINGLY. FORMAT: 



CALL PLOTD(XDATA, YDATA, NNPTS, EJECT, LABEL, XLABEL, YLABEL, = 
X MINMAX, CRVTTL, CRVNUM, XMIN, XMAX, YMIN, YMAX) 

WHERE: 

XDATA IS A REALXA ARRAY DIMENSIONED AT LEAST | NNPTS | 

CONTAINING THE X ORDINATE VALUES, 

YDATA IS A REALXA ARRAY DIMENSIONED AT LEAST | NNPTS | = 

CONTAINING THE Y ORDINATE VALUES, 

NNPTS IS AN INTEGERXA SCALAR DESIGNATING THE NUMBER OF = 
POINTS TO BE PLOTTED. THE NUMBER OF POINTS IS 
ABS(NNPTS). NNPTS<0 MEANS PLOT POINTS ONLY. 

EJECT IS A LOGICALXA VARIABLE OR CONSTANT INDICATING = 
WHETHER A PAGE EJECT IS REQUIRED FOLLOWING THE = 
CURRENT CURVE. THIS ALLOWS MULTIPLE CURVES ON ONE = 
SET OF EXES. PAGE EJECT WILL OCCUR FOR NEXT GRAPH = 
AFTER EJECT HAS BEEN SET TO .TRUE. 



LABEL IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY = 

CONTAINING THE INTENDED LABEL FOR THE GRAPH. THE = 
MAXIMUM ALLOWABLE LENGTH (INCLUDING CHARACTER) = 
IS 32 CHARACTERS. 

XLABEL IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY = 

CONTAINING THE INTENDED LABEL OF THE X-AXIS OF THE = 
GRAPH. IN THIS PROGRAM, XLABEL IS ALWAYS 'ALPHA'. = 

YLABEL IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY = 

CONTAINING THE INTENDED LABEL OF THE Y-AXIS OF THE = 
GRAPH. IN THIS PROGRAM, YLABEL IS ALWAYS 'BETA*. 

MINMAX IS A PARAMETER THAT DETERMINES WHETHER THE MINIMUM = 
AND MAXIMUM VALUES FOR THE AXES ARE TO BE ASSIGNED = 
BY THE USER, OR WHETHER THEY WILL BE 'AUTOSCALED' . = 

CRVTTL IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY = 

AND TERMINATED BY A CHARACTER SPECIFYING THE = 

INTENDED NAME WHICH LABELS AN INDIVIDUAL CURVE. 

CRVNUM IS A REAL VARIABLE OR CONSTANT THAT SPECIFIES THE = 
VALUE TO BE CONCATENATED ONTO THE END OF 'CRVTTL'. = 
FOR EXAMPLE, IF THIS CURVE REPRESENTS 'ZETA =0.5' = 

THEN CRVTTL = 'Z=$', WHILE CRVNUM = 0.5. 



C 



SUBROUTINE PLOTD(XDATA, YDATA, NNPTS, EJECT, LABEL, 
X YLABEL, MINMAX, CRVTTL, CRVNUM, 

X XMIN, XMAX, YMIN, YMAX,GRD) 

REALXA XDATA(l), YDATA(l) 

REALXA XMIN, XMAX, YMIN, YMAX 

REALXA CRVTTL, CRVNUM 

INTEGERXA NPTS, NNPTS 

LOGICALXA EJECT 

LOGICALXI LABEL(l) 

LOGICALXl XLABEL(l) 

L0GICAL5(1 YLABEL(l) 



XLABEL, 



LOGICALXA INIT /.FALSE./ 
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SET THE NUMBER OF POINTS 
NPTS = lABS(NNPTS) 

IF THE ROUTINE HAS BEEN NOT BEEN INITIALIZED (INIT = .FALSE.) 
THEN INITIALIZE IT. 

IF (.NOT. INIT) CALL PLDOOl (XDATA, YDATA, NPTS, LABEL, XLABEL, 

X YLABEL, MINMAX, CRVTTL , CRVNUM, 

X XMIN, XMAX, YMIN, YMAX,GRD) 

FRAME THE PLOT AND REORDER SYMBOLS 

IF (.NOT. INIT) CALL FRAME 

INDICATE INITIALIZATION IN CASE THIS ROUTINE IS RE-ENTERED 
WITH MULTIPLE CURVES. 

INIT = .TRUE. 

CALCULATE THE NUMBER OF POINTS TO GIVE ABOUT 5 MARKERS PER 
LINE 

TENTATIVELY SET POINTS ONLY, THEN CHECK 
NMARK = -1 

IF(NNPTS.LE.O)GO TO 10 

CURVE WANTED, SET MARKERS 
IF (MOD(NPTS, A) .EQ. 1) NMARK = NPTS / A 
IF (MOD(NPTS, A) .NE. 1) NMARK = NPTS / 3 
IF (NMARK .EQ. 0) NMARK = 1 
10 CONTINUE 



DRAW THE CURVE 

XXMAX = -1.0E75 
YYMAX = -1.0E75 

DO 20 1=1, NPTS 

IF(XDATAd) .GT. XMAX .OR. XDATA(I) . LT . XMIN .OR. 

X YDATA(I) .GT. YMAX .OR. YDATA(I) .LT. YMIN) GO TO 20 

IF(XXMAX .LT. XDATA(D) J = I 

IF(XXMAX .LT. XDATA(D) XXMAX = XDATA(I) 

IF(YYMAX .LT. YDATA(D) K=I 

IF(YYMAX .LT. YDATA(D) YYMAX = YDATA(I) 

20 CONTINUE 

IF((YMAX-YYMAX) - (XMAX-XXMAX) ) AO, 30, 30 
30 L=J 

GO TO 50 
AO L=K 

50 CALL GRACE(0.) 

CALL DOT 

IF (GRD) 60,60,70 
60 CALL GRID(1,1) 

70 CALL RESET! ’DASH* ) 

IF (CRVNUM .EQ. -9.7531) GO TO 80 
CALL HEIGHT(0.125) 

CALL RLMESS(CRVTTL , 3,XDATA( L ) ,YDATA( L ) ) 

CALL RLREAL(CRVNUM,2, ’ABUT’, ’ABUT’) 

CALL THKCRV(0.015) 

80 CALL CURVE(XDATA, YDATA, NPTS, NMARK) 

CALL RESET! ’THKCRV’) 

IF THIS IS NOT THE LAST (OR ONLY) CURVE ON THIS GRAPH, THEN 
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EXIT. OTHERWISE CLOSE THE PLOT AND TURN OFF INITIALIZATION 
FLAG. 

IF (.NOT. EJECT) RETURN 

END OF THIS PLOT 

CALL ENDPL(O) 

INIT = .FALSE. 

RETURN 

END 



SUBROUTINE PLDOOKXDATA, YDATA, NPTS, LABEL, XLABEL, 
X YLABEL, MINMAX, CRVTTL, CRVNUM, 

X XMIN, XMAX, YMIN, YMAX, GRD) 

THIS SUBROUTINE DOES THE INITIALIZING. 



REAL3^4 


XDATA(NPTS) 


REALXA 


YDATA(NPTS) 


INTEGER5(^ 


NPTS 


LOGICALXl 


LABEL(l) 


LOGICALXl 


XLABEL(l) 


LOGICALXl 


YLABEL(l) 


REAL5(A 


XMIN 


REALX4 


XMAX 


REAL5(A 


YMIN 


REAL^^*^ 


YMAX 



INITIALIZE DISSPLA 
CALL PLD009 

CALL HEADIN(LABEL, 100, 2., 1) 
CALL XNAMECXLABEL, 100) 

CALL YNAMECYLABEL, 100) 



EXTRACT MINIMA AND MAXIMA 

IF (MINMAX .NE. 1) GO TO 90 

CALL PLD010(XDATA, NPTS, XMIN, XMAX) 

CALL PLD010(YDATA, NPTS, YMIN, YMAX) 

CALL THE LINEAR-LINEAR INITIALIZING ROUTINE 

90 CALL PLDOll (XMIN, XMAX, YMIN, YMAX) 

RETURN 

END 



SUBROUTINE PLD009 

THIS SUBROUTINE ESTABLISHES THE PARAMETERS FOR DISSPLA. 

NOTE THAT IT IS THE USER'S RESPONSIBILITY TO NOMINATE THE 
GRAPHIC DEVICE. 



CALL NOCHEK 
CALL NOBRDR 
CALL PAGEdA. ,1A. ) 
CALL PHYS0R(2., .75) 

GO TO LEVEL 2. 

CALL AREA2D(9.,9. ) 



107 



O O o ooooo ooo ooo oo ooo ooo ooooooooooooo oo ooo oo 



LETTERING IS DUPLX WITH UPPER CASE ONLY. 

CALL DUPLX 

CALL BASALF( 'STANDARD') 

INTEGER (OR ROUNDED) AXES WITH Y AXIS LABELLING AT 0 DEGREES. 

CALL INTAXS 
CALL YAXANG(0.) 

RETURN 

END 



SUBROUTINE PLD0I0(V, N, MIN, MAX) 

THIS SUBROUTINE SCANS A VECTOR FOR MAXIMUM AND MINIMUM 
INPUT PARAMETERS: 

V DATA VECTOR (REAL) 

N NUMBER OF POINTS IN VECTOR V (INTEGER) 

OUTPUT PARAMETERS: 



MIN VECTOR ORIGIN (REAL) 

MAX VECTOR MAXIMUM (REAL) 



REALX^i 


V(N) 


INTEGERX4 


N 


REALXA 


MIN 


REALX4 


MAX 



INITIALIZE THE MAXIMA AND MINIMA 



MIN = I.0E75 
MAX = -1.0E75 

FIND MAXIMUM AND MINIMUM OF VECTOR V 



DO 100 I = 1, N 
IF (MIN .GT. V(D) MIN = V(I) 
IF (MAX .LT. V(D) MAX = V(I) 
100 CONTINUE 
RETURN 
END 



SUBROUTINE PLDOll (XMIN,XMAX, YMIN, YMAX) 

THIS SUBROUTINE SETS UP DISSPLA FOR A LINEAR-LINEAR AXIS PLOT. 

REALXA XMIN 

REALX^i XMAX 

REALXA YMIN 

REALXA YMAX 

A SIMPLE CALL TO GRAF WILL DO IT... 

CALL HEIGHT(0.175) 

CALL GRAF(XMIN, 'SCALE', XMAX, YMIN, 'SCALE', YMAX) 

RETURN 

END 
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SUBROUTINE ROOTS — CALCULATES ROOTS OF THE NO ORDER EQUATION 
FOR EVERY TENTH VALUE ALPHA/BETA PAIR. 



SUBROUTINE ROOTS (ALPHA, BETA, AJ, BJ, CJ, N02) 
INTEGERBCA N02, NN, NNN, NNNN 

REAL5(A ALPHA, BETA, AJ(IOO), BJ(IOO), CJ(IOO) 

REAL^^S COEF(IOO), ROOTMY(IOO) 

WRITE(6, AO) 

NN = N02+1 
NNNN = N02 
DO 10 1=1, NN 
NNN = NN+l-I 

10 COEF(NNN) = CJ(I) + (ALPHA AJ(D) + (BETA ^ BJ(D) 
CALL ZPOLR (C0EF,N02,R00TMY,IER) 

WRITE (6,50) 

DO 20 1=1, NNNN 

II = 2X1 

III = II-l 

20 WRITE (6,60) ROOTMY(III), ROOTMY(II) 

WRITE (6,30) 

30 FORMAT(//////////////////// ) 

AO FORMAT(/20X, * ROOTS FOR ABOVE ALPHA, BETA’) 

50 F0RMAT(21X,* REAL PART IMAG PART’) 

6 0 FORMAT (21X,E10.A,6X,E10.A) 

RETURN 

END 



SUBROUTINE SAVIT -- SAVES DATA IN FN FT FM = INAME DATA Al, 
WHERE INAME IS THE USER’S CHOICE. 



SUBROUTINE SAVIT 

COMMON /SAVE/ LABEL, WN, ND, N02, NC, C J , AJ, BJ, 

X NZ, ZETA, NS, SIGMA, NW, W, NZW, ZW, 

X XMIN, XMAX, YMIN, YMAX 

REALXA WN, CJ(IOO), AJ(IOO), BJ(IOO), XMIN, XMAX, YMIN, YMAX 
REALXA ZETA(IOO), SIGMA(IOO), W(IOO), ZW(IOO) 

INTEGER ND, N02, NC, NZ, NS, NW, NZW 
CHARACTERXA LABEL(9), INAME(2) 

WRITE(6,10) 



READ(5,20) (INAME(I), 1=1,2) 
CALL FRTCMSt 'FILEDEF ','02 




', 'DISK 


', INAME, 'DATA 


') 


WRITE(2,30) (LABEL(I), 1=1, 
WR1TE(2,X) NN 

HR1TE(2,X) ND, N02, NC, NZ, 


9) 

NS, 


NN, NZN 







WRITE(2,x) (CJ(J), J=l, NC) 

WRITE(2,X) (AJ(J), J=l, NC) 

WRITE(2,X) (BJ(J), J=l, NC) 

WRITE(2,x) (ZETA(M), M=l, NZ) 

WRITE(2,x) (SIGMA(M), M=l, NS) 

WRITE(2,x) (W(M), M=l, NW) 

WRITE(2,x) (ZW(M), M=l, NZW) 

WRITE(2,X) XMIN, XMAX, YMIN, YMAX 
10 F0RMAT(5X, ’UNDER WHAT NAME DO YOU WANT TO SAVE THE DATA?’,/, 
X 5X,’(8 CHARACTERS MAX)’) 

20 F0RMAT(2AA) 

30 F0RMAT(9A4) 

END FILE 02 
REWIND 02 
RETURN 
END 
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SUBROUTINE GETIT — RETRIEVES DATA FROM FN FT FM = INAME DATA Al, = 
WHERE INAME IS THE USER'S CHOICE. 



SUBROUTINE GETIT 

COMMON /SAVE/ LABEL, WN, ND, N02, NC, CJ, AJ, BJ, 

X NZ, ZETA, NS, SIGMA, NW, W, NZW, ZW, 

X XMIN, XMAX, YMIN, YMAX 

REALXA WN, CJ(IOO), AJ(IOO), BJ(IOO), XMIN, XMAX, YMIN, YMAX 
REALXA ZETA(IOO), SIGMA(IOO), W(IOO), ZW(IOO) 

INTEGER ND, N02, NC, NZ, NS, NW, NZW 
CHARACTERX4 LABEL(9) 

CHARACTERX8 INAME 
CHARACTERX21 NAME 
10 WRITE(6,40) 

READ(5,50) INAME 

NAME = 'STATE '//INAME//' DATA X' 

CALL EXCMS(NAME,RC) 

IF (RC .EQ. 0) GOTO 20 
WRITE(6,30) 

GOTO 10 

20 CALL FRTCMS( 'FILEDEF ','02 ','DISK ', INAME, ' DATA ') 

READ(2,60) (LABEL(I), 1=1, 9) 

READ(2,X) WN 

READ(2,X) ND, N02, NC, NZ, NS, NW, NZW 
READ(2,X) (CJ(J), J=l, NC) 

READ(2,X) (AJ(J), J=l, NC) 

READ(2,X) (BJ(J), J=l, NC) 

READ(2,X) (ZETA(M), M=l, NZ) 

READ(2,X) (SIGMA(M), M=l, NS) 

READ(2,X) (W(M), M=l, NW) 

READ(2,X) (ZW(M), M=l, NZW) 

READ(2,X) XMIN, XMAX, YMIN, YMAX 
30 F0RMAT(5X, 'DATA FILE NOT FOUND - TRY ANOTHER') 

40 F0RMATC5X, 'UNDER WHAT NAME IS THE DATAFILE SAVED? ',/, 

X 5X, '(8 CHARACTERS MAX)') 

50 FORMATdAS) 

60 F0RMAT(9A4) 

END FILE 02 
REWIND 02 
RETURN 
END 



SUBROUTINE ASTER PLACES A $ AT THE END OF A CHARACTER STRING 



SUBROUTINE ASTERCLLINES, LINE) 

CHARACTERX4 DOLLAR/'^ '/, LLINES(8), LINE(9) 
DO 10 1=1,8 

LINE(I)=LLINES(I) 

CONTINUE 

LINE(9)=D0LLAR 

RETURN 

END 
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c 

c ==================================================================== 

C SUBROUTINE READR — INTERACTIVELY READS A REAL NUMBER REPLY. 

C IF THE USER INADVERTENTLY ENTERS A NULL STRING 
C A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. 

C ==================================================================== 

C 

C 

C 

SUBROUTINE READR (ANSR) 

REALXA ANSR 
INTEGER COUNT 
COUNT=0 

10 CONTINUE 

COUNT=COUNT+l 
IF (COUNT.lt. 3) GO TO 20 
WRITE (5,60) ■ 

GO TO AO 
20 CONTINUE 

READ (5,x,END=30,ERR=30) ANSR 
RETURN 

30 REWIND 5 

WRITE (5,50) 

GO TO 10 
AO CONTINUE 

STOP 

50 FORMAT (IX, » WARNING: NULL STRINGS ARE NOT ALLOWED, ENTER A NUMER 

XICAL VALUE. ') 

60 FORMAT (///,5X,* PROGRAM TERMINATION - TWO NULL STRINGS ENTERED!’) 

END 
C 
C 

C ==================================================================== 

C SUBROUTINE READI -- INTERACTIVELY READS AN INTEGER REPLY. 

C IF THE USER INADVERTENTLY ENTERS A NULL STRING OR NEGATIVE VALUE = 

C A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. 

C ==================================================================== 

C 

C 

SUBROUTINE READI (IANS) 

INTEGER COUNT, IANS 
COUNT=0 

10 CONTINUE 

COUNT=COUNT+l 
IF (COUNT.lt. 3) GO TO 20 
WRITE (5,70) 

GO TO 50 
20 CONTINUE 

READ (5,x,END=A0,ERR=A0) IANS 
IF (IANS) 40,30,30 
30 CONTINUE 

RETURN 

40 REWIND 5 

WRITE (5,60) 

GO TO 10 
50 CONTINUE 

STOP 

60 FORMAT (IX,’ WARNING: IMPROPER DATA ENTRY! ENTER A POSITIVE INTEG 
XER. ’ ) 

70 FORMAT (///,5X,’ PROGRAM TERMINATION - 2 IMPROPER DATA ENTRIES!’) 

END 
C 
C 
C 
C 
C 
C 
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c 

c ==================================================================== 

C SUBROUTINE READC " INTERACTIVELY READS A CHAR STRING REPLY. 

C (»YES» OR *NO»). IF THE USER INADVERTENTLY ENTERS A NULL STRING 
C A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. 

C ==================================================================== 

C 

C 

SUBROUTINE READC (CANS) 

INTEGER COUNT 
CHARACTERX4 CANS 
COUNT=0 

10 CONTINUE 

COUNT=COUNT+l 
IF (COUNT.lt. 3) GO TO 20 
WRITE (5.60) 

GO ‘TO AO 
20 CONTINUE 

REWIND 5 

READ (5.70.END=30.ERR=30) CANS 
RETURN 

30 REWIND 5 

WRITE (5.50) 

GO TO 10 
40 CONTINUE 

STOP 

50 FORMAT (IX. » WARNING: NULL STRINGS ARE NOT ALLOWED ») 

60 FORMAT (///.5X.» PROGRAM TERMINATION - TWO NULL STRINGS ENTERED! ») 

70 FORMAT (A2) 

END 

C 

C 

C ==================================================================== 

C SUBROUTINE READL — INTERACTIVELY READS A STRING OF CHARACTERS. 

C IF THE USER INADVERTENTLY ENTERS A NULL STRING 
C A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. 

C ==================================================================== 

C 

C 

SUBROUTINE READL( LLINES) 

INTEGER COUNT. I. NIX 

CHARACTER^A BBLANK/» V, LLINES(8) 

DO 10 1=1.8 

LLINES(I) = BBLANK 
10 CONTINUE 
C0UNT=0 

20 COUNT=COUNT+l 

IF(C0UNT.LT.3) GO TO 30 
WRITE(6.70) 

GO TO 50 
30 CONTINUE 
REWIND 5 

READ(5.80. END=A0. ERR=30)(LLINES(J). J=1.9) 

RETURN 

AO REWIND 5 

WRITE(6.60) 

GO TO 20 
50 CONTINUE 
STOP 

60 F0RMAT(1X.» WARNING: NULL STRINGS ARE NOT ALLOWED. ENTER CHARACTER 
X VALUES. ») 

70 FORMAT (///.5X.» PROGRAM TERMINATION - TWO NULL STRINGS ENTERED! ») 
80 F0RMAT(9AA) 

END 

C 

C 

C 

C 

C 

C 



112 



oooooooo 



SUBROUTINE META — BY CALLING COMPRS AS A SUBROUTINE HERE. 
DONEPL HAS SUFFICIENT TIME TO FINISH; OTHERWISE COMPRS IGNORED 
AND METAFILE NOT SAVED 



SUBROUTINE META 
DO 10 1=1,900000 
10 CONTINUE 

CALL COMPRS 

RETURN 

END 
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